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ABSTRACT 

This  paper  has  two  main  objectives.  First,  it  develops  the  multifractal  formalism  in  a 
context  suitable  for  both,  measures  and  functions,  deterministic  as  well  as  random,  thereby 
emphasizing  an  intuitive  approach.  Second,  it  carefully  discusses  several  examples,  such  as 
the  binomial  cascades  and  self-similar  processes  with  a  special  eye  on  the  use  of  wavelets. 
Particular  attention  is  given  to  a  novel  class  of  multifractal  processes  which  combine  the 
attractive  features  of  cascades  and  self-similar  processes.  Statistical  properties  of  estimators 
as  well  as  modelling  issues  are  addressed. 
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1  Introduction  and  Summary 

Fractal  processes  have  been  instrumental  in  a  variety  of  fields  ranging  from  the  theory  of 
fully  developed  turbulence  [73,  64,  36,  12,  7],  to  stock  market  modelling  [28,  68,  69,  80], 
image  processing  [61,  21,  104],  medical  data  [2,  98,  11]  and  geophysics  [36,  65,  47,  92],  In 
networking,  models  using  fractional  Brownian  motion  (fBm)  have  helped  advance  the 
held  through  their  ability  to  assess  the  impact  of  fractal  features  such  as  statistical  self¬ 
similarity  and  long-range  dependence  (LRD)  to  performance  [60,  81,  90,  89,  96,  34,  88]. 

Roughly  speaking,  a  fractal  entity  is  characterized  by  the  inherent,  ubiquitous  oc¬ 
currence  of  irregularities  which  governs  its  shape  and  complexity.  The  most  prominent 
example  is  certainly  fBm  Bu(t)  [71].  Its  paths  are  almost  surely  continuous  but  not 
differentiable.  Indeed,  the  oscillation  of  fBm  in  any  interval  of  size  5  is  of  the  order  5H 
where  H  e  (0, 1)  is  the  self-similarity  parameter: 

BH(at)  =  aHBH(t).  (1.1) 

Reasons  for  the  success  of  fBm  as  a  model  of  LRD  may  be  seen  in  the  simplicity  of  its 
scaling  properties  which  makes  it  amendable  to  analysis.  The  fact  of  being  Gaussian 
bears  further  advantages.  However,  the  scaling  law  (1.1)  implies  also  that  the  oscil¬ 
lations  of  fBm  at  fine  scales  are  uniform*  which  comes  as  a  disadvantage  in  various 
situations  (see  Figure  1).  Real  world  signals  often  possess  an  erratically  changing  os¬ 
cillatory  behavior  (see  Figure  2)  which  have  earned  them  the  name  multifractals,  but 
which  also  limits  the  appropriateness  of  fBm  as  a  model.  This  rich  structure  at  fine 
scales  may  serve  as  a  valuable  indicator,  and  ignoring  it  might  mean  to  miss  out  on 
relevant  information  (see  references  above). 

This  paper  has  two  objects.  First,  we  present  the  framework  for  describing  and 
detecting  such  a  multifractal  scaling  structure.  Doing  so  we  survey  local  and  global 
multifractal  analysis  and  relate  them  via  the  multifractal  formalism  in  a  stochastic 
setting.  Thereby,  the  importance  of  higher  order  statistics  will  become  evident.  It  might 
be  especially  appealing  to  the  reader  to  see  wavelets  put  to  efficient  use.  We  focus 
mainly  on  the  analytical  computation  of  the  so-called  multifractal  spectra  and  on  their 
mutual  relations.  Thereby,  we  emphasize  issues  of  observability  by  striving  for  formulae 
which  hold  for  all  or  almost  all  paths  and  by  pointing  out  the  necessity  of  oversampling 
needed  to  capture  certain  rare  events.  Statistical  properties  of  estimators  of  multifractal 
quantities  as  well  as  modelling  issues  are  addressed  elsewhere  (see  [41,  3,  40]  and 
[68,  89,  88]). 

Second,  we  carefully  discuss  basic  examples  as  well  as  Brownian  motion  in  multi¬ 
fractal  time ,  Bi/2(M.{t)).  This  process  has  recently  been  suggested  as  a  model  for  stock 
market  exchange  by  Mandelbrot  who  argues  that  oscillations  in  exchange  rates  occur 
in  multifractal  ‘trading  time’  [68,  69].  With  the  theory  developed  in  this  paper,  it  be¬ 
comes  an  easy  task  to  explore  Bi/2(M.{t))  from  the  multifractal  point  of  view,  and  with 


*This  property  is  also  known  as  the  Lvey  modulus  of  continuity  in  the  case  of  Brownian 
motion.  For  fBm  see  [5,  Thm.  8.3.1.]. 
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FIGURE  1.  Fractional  Brownian  motion,  as  well  as  its  increment  process  called  fGn  (displayed 
on  top  in  T5),  has  only  one  singularity  exponent  h(t)  =  H,  a  fact  which  is  represented  in 
the  linear  partition  function  r  (see  T2)  and  a  multifractal  spectrum  (see  T3)  which  consists 
of  only  one  point:  for  fBm  (H,  1)  and  for  fGn  ( H  —  1,1).  For  further  details  on  the  plots  see 
(1.9),  (1.6)  and  Figure  7. 

little  more  effort  also  more  general  multifractal  ‘subordinators’.  The  reader  interested 
in  these  multifractal  processes  may  wish,  at  least  at  first  reading,  to  content  himself 
with  the  notation  introduced  on  the  following  few  pages,  skip  the  sections  which  deal 
more  carefully  with  the  tools  of  multifractal  analysis,  and  proceed  directly  to  the  last 
sections.  The  remainder  of  this  introduction  provides  a  summary  of  the  contents  of  the 
paper,  following  roughly  its  structure. 

1.1  Singularity  Exponents 

In  this  work,  we  are  mainly  interested  in  the  geometry  or  local  scaling  properties 
of  the  paths  of  a  process  Y(t).  Therefore,  all  notions  and  results  concerning  paths 
will  apply  to  functions  as  well.  The  study  of  fine  scale  properties  of  functions  (as 
opposed  to  measures)  has  been  pioneered  in  the  work  of  Arneodo,  Baery  and  Muzy 
[7,  78,  79,  1,  2,  80],  who  were  also  the  first  to  introduce  wavelet  techniques  in  this 
context.  For  simplicity  of  the  presentation  we  take  t  G  [0, 1].  Extensions  to  the  real  line 
1R  as  well  as  to  higher  dimensions,  being  straightforward  in  most  cases,  are  indicated. 

A  typical  feature  of  a  fractal  process  Y(t)  is  that  it  has  a  non-integer  degree  of 
differentiability,  giving  rise  to  an  interesting  analysis  of  its  local  Holder  exponent  H(t) 
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FIGURE  2.  Real  world  signals  such  as  this  geophysical  data  often  exhibit  erratic  behavior 
and  their  appearance  may  make  stationarity  questionable.  One  such  feature  are  ‘trends’ 
which  sometimes  can  be  explained  by  strong  correlations  (LRD).  Another  such  feature  are 
the  sudden  jumps  or  ‘ bursts’  which  in  turn  are  a  typical  for  multifractals.  For  such  signals 
the  singularity  exponent  h(t )  depends  erratically  on  time  t,  a  fact  which  is  reflected  in  the 
concave  partition  function  r  (see  T2)  and  a  multifractal  spectrum  (see  T3)  which  extends 
over  a  non-trivial  range  of  singularity  exponents. 

which  is  roughly  defined  through 

\Y(t')  -P(t')\  ~  \t'^t\H{t)  (1.2) 

for  some  polynomial  P  which  in  nice  cases  is  simply  the  Taylor  polynomial  of  Y  at  t. 
A  rigorous  definition  is  given  in  (2.1). 

Provided  the  polynomial  is  constant,  H(t)  can  be  obtained  from  the  limiting  behavior 
of  the  so-called  coarse  Holder  exponents,  i.e., 

he(t)  =  y-— log  sup  \Y(t')  -  Y(t)\.  (1.3) 

logs  |t'— 1|<£ 

For  rigorous  statements  we  refer  to  (2.2)  and  lemma  2.3. 

However,  as  the  example  t2+t2'‘  shows,  the  use  of  h£(t)  is  ineffective  in  the  presence  of 
polynomial  trends.  Then,  he(t)  will  reflect  the  lowest  non-constant  term  of  the  Taylor 
polynomial  of  Y  at  t.  For  this  reason,  and  also  to  avoid  complications  introduced 
through  the  computation  of  the  supremum  in  (1.3),  one  may  choose  to  employ  wavelet 
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decompositions  or  other  tools  of  time  frequency  analysis.  Properly  chosen  wavelets  are 
blind  to  polynomials  and  due  to  their  scaling  properties  they  contain  information  on 
the  Holder  regularity  of  Y  [51,  23].  Their  application  in  multifractal  estimation  has 
been  pioneered  by  [7,  53,  30].  Furthermore,  wavelets  provide  unconditional  basis  for 
several  regularity  spaces  such  as  Besov  spaces  (see  (2.14)  and  (6.2))  whence  their  use 
bears  further  advantages. 

Yet,  the  ‘classical’  choice  of  a  singularity  exponent  is 

<4”’  =  (M((k+  1)2-”)  -M(k2-")).  (1.4) 

It  is  attractive  due  to  its  simplicity  and  becomes  actually  quite  powerful  when  studying 
monotonously  increasing  processes  M(t),  in  particular  the  distribution  functions  of 
singular  measures,  such  as  cascades. 

In  this  chapter  we  will  introduce  the  exponents  w ^  emerging  from  a  wavelet  based 
analysis  and  elaborate  on  the  relation  between  these  different  singularity  exponents 
h^\  and  w^\ 

1.2  Multifractal  Spectra 

As  indicated  we  are  mainly  interested  in  the  geometry  or  local  regularity  of  the  paths 
of  Y{t).  Let  us  fix  such  a  realization  for  the  time-being. 

Local  analysis 

Ideally,  one  would  like  to  quantify  in  geometrical  as  well  as  statistical  sense  which 
values  H(t )  appear  on  a  given  path  of  the  process  Y,  and  how  often  one  will  encounter 
them.  Towards  the  first  description  one  studies  the  sets 

4“'  =  {*  :  H(t)  =  a}  (1.5) 

for  varying  a.  Similarly,  one  could  consider  sets  E$  and  E$  defined  through  the 
limiting  behavior  of  the  singularity  exponents  cxf'>  or  respectively.  If  no  confusion 
regarding  the  choice  of  /qn) ,  u/f1'1  or  a ^  can  arise,  we  simply  drop  the  index. 

The  sets  E^a  form  a  decomposition  of  the  support  of  Y  according  to  its  singularity 
exponents.  We  say  that  Y  has  a  rich  multifractal  structure  if  these  sets  E^  are  highly 
interwoven,  each  lying  dense  on  the  line.  Typically,  only  one  of  the  E^  has  full  Lebesgue 
measure,  while  the  others  form  dusts,  more  precisely,  sets  with  non-integer  Hausdorff 
dimension  dim(i7[al)  [32],  Dimensions  are  always  positive,  and  the  smaller  the  dimension 
of  a  set  the  ‘thinner’  the  set.  In  this  sense,  the  function 

a  i — *  dim(FJ^)  (1.6) 

gives  a  compact  representation  of  the  complex  singularity  structure  of  Y .  It  has  been 
termed  the  multifractal  spectrum  of  Y  and  is  studied  extensively  in  the  ‘classical’  liter¬ 
ature. 
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To  develop  some  intuition  let  us  consider  a  differentiable  path.  To  avoid  trivialities 
let  us  assume  that  this  path  and  its  derivative  have  no  zeros.  Then,  dim ( A  A  )-spectrum 
reduces  to  the  point  (1, 1).  On  the  other  hand,  if  H(t)  is  continuous  and  not  constant 
on  intervals  then  each  A^  is  finite  and  dim(A[“])  =  0  for  all  a  in  the  range  of  H{t). 
A  spectrum  dirn(AA)  with  non-degenerate  form  is,  thus,  indeed  indication  for  rich 
singularity  behavior.  By  this  we  mean  that  H{t)  changes  erratically  with  t  and  takes 
each  value  a  on  a  rather  large  set  A  A. 


Global  analysis 

A  simpler  way  of  capturing  the  complex  structure  of  a  signal  is  obtained  when  adapting 
the  concept  of  box-dimension  to  the  multifractal  context.  As  the  name  indicates,  one 
aims  at  an  estimate  of  dim(i?[al)  by  counting  the  intervals  -  or  boxes  -  over  which  Y 
increases  roughly  with  the  ‘right’  Holder  exponent.  Therefore,  we  need  to  introduce 
grain  exponents,  a  discrete  approximation  to  h£{t)  (see  (1.3)): 

—  (l/n)log2sup{|y(s)  —  F(t)|  :  (/c  —  l)2_n  <  s  <  1  <  (/c  +  2)2“n}  (1.7) 


and  define  the  grain  (multifractal)  spectrum  as  [73,  46,  45,  91] 

logAfW(a,  e) 


/(«)  —  iimlimsup 

£  ,0  77, — »oO  ^  2 


(1.8) 


where  N^n\a,£)  =  #{k  :  | —  a|  <  e}  counts,  how  many  of  the  grain  exponents  h ^ 
are  approximately  equal  to  a.  Similarly,  one  may  define  such  spectra  for  the  singularity 
exponents  a ^  and  If  confusion  may  arise,  we  will  indicate  the  chosen  exponent 

by  writing  explicitly  fh(a),  /Q(a),  or  /„,(«). 

This  multifractal  spectrum  can  be  interpreted  (at  least)  in  three  ways.  First,  as 
mentioned  already  it  is  related  to  the  notion  of  dimensions.  Indeed,  a  simple  argument 
shows  that  dim(A[a])  <  /(a)  [94],  The  essential  ingredient  for  a  proof  is  the  fact  that  the 
calculation  of  dim  (A  A)  involves  finding  an  optimal  covering  of  A  A  while  /(a)  considers 
only  uniform  covers.  In  short,  /(a)  provides  an  upper  bound  on  the  dimension  and, 
thus,  the  ‘size’  of  the  sets  A  A . 

Second,  (1.8)  suggests  that  the  re-normalized  histograms  (1/n)  log2  N^n\a,  e)  should 
all  be  roughly  equal  at  small  scales  2~n  to  the  scale  independent  /(a).  It  should  be 
remembered  that  this  is  foremost  (by  definition)  a  property  of  the  paths  of  the  given 
process.  We  stress  this  point  because  it  is  tempting  to  argue  that  -at  least  under  suitable 
ergodicity  assumptions-  one  should  see  the  marginal  distribution  of  }if']  reflected  in  /. 
However,  one  should  not  overlook  that  the  logarithmic  re-normalization  implemented  in 
/(a)  is  aimed  at  detecting  exponential  scaling  properties  rather  than  the  marginals  on 
multiple  scales  themselves.  For  fBm  (see  (1.1))  this  re-normalization  indeed  causes  all 
details  of  the  Normal  multi-scale  marginals  to  be  washed  out  into  a  virtually  structure¬ 
less  /(a)  which  gives  notice  of  the  presence  of  only  one  scaling  law,  the  self-similarity 
(1.1)  with  parameter  H.  Thus,  /  expresses  that  fBm  is  ‘mono-fractal’,  as  mentioned 
above.  To  the  contrary  with  ‘multi-fractal’  processes  such  as  multiplicative  cascades, 
for  which  /  reflects  the  presence  of  an  entire  range  of  scaling  exponents  (see  (5.32)). 


1  Introduction  and  Summary  7 


The  third  natural  context  for  the  coarse  spectrum  /  is  that  of  Large  Deviation 
Principles  (LDP)  [29,  91].  Indeed,  N^n\a,£)/ 2”  defines  a  probability  distribution f  on 
{ tif'1  :  k  =  0, ...  ,2”  —  1}.  Alluding  to  the  Law  of  Large  Numbers  (LLN)  we  may 
expect  this  distribution  to  be  concentrated  more  and  more  around  the  ‘most  typical’ 
or  ‘expected’  value  as  n  increases.  The  spectrum  /(a)  measures  how  fast  the  chance 
N('n\a,£)/2n  to  observe  a  ‘deviant’  value  a  decreases,  i.e.,  N^n\a,  e)/2n  ~  2-Ra)_1. 

The  close  connection  to  LDP  leads  one  to  study  the  scaling  of  ‘sample  moments’ 
through  the  so-called  partition  function  [45,  46,  36,  91] 

1  c(n)/  x  2"  — 1 

rfe(g)  :=  liminf  S  where  S^\q)  :=  V"  2~nqh<?\  (1.9) 

n—>oo  —n  log  2  z — ' 

k= 0 

which  are  defined  for  all  q  E  IR.  Similarly,  replacing  hY  by  ,  one  defines  ra{q)  and 
Sa(q).  The  latter  takes  on  the  well-known  form  of  a  partition  sum 

2n  — 1 

Sj?\q)  =  2~nqa*)  =  \Y  {(k  +  1)2~n)  ~  Y{k2~n)\q  .  (1.10) 

k= 0 

Again  similarly,  one  defines  Tw{q)  and  Sw\q )  by  replacing  Iil[l>  by  wavelet  based  ex¬ 
ponents  w (see  (2.11)).  Again,  if  no  confusion  on  the  choice  of  h\f\  wY  or  aY  can 
arise,  we  simply  drop  the  index  h,  a  or  w  . 


1.3  Multifractal  Formalism 

The  theory  of  LDP  suggests  /(a)  and  r{q)  are  strongly  related  since  2~nS<'n\q )  is  the 
moment  generating  function  of  the  random  variable  An{k )  :=  *  ln(2)  (recall  foot¬ 

note  f).  For  a  motivation  of  a  formula  connection  /(a)  and  r(g)  consider  the  heuristics 


S(n\q )  = 


EE 


.  (n) 

2  ~nihk  ~ 


2nfA)2~n(ia  = 


2~n(qa-f(a))  ^  2~n  inf“(9“-/(“)) 


j  (n) 

K  -a 


Assuming  that  'Yha  has  only  finite  many  terms  the  last  step  simply  replaces  the  sum  by 
its  strongest  term.  Making  this  entire  argument  rigorous  we  prove  in  this  paper  that 


T(q)  =  f*(a )  :=  mi(qa  -  /(a)).  (1.11) 

a 

Here  (•)*  denotes  the  Legendre  transform  which  is  omnipresent  in  the  theory  of  LDP. 
Indeed,  by  applying  a  theorem  due  to  Gartner  and  Ellis  [27]  and  imposing  some  reg¬ 
ularity  on  r(g)  theorem  3.5  shows  that  the  family  of  probability  densities  defined  by 
N(n\a,e)/ 2n  satisfies  the  full  LDP  [26]  with  rate  function  f  meaning  that  /  is  actually 
a  double- limit  and  /(a)  =  r*(a).  Corollary  4.5  establishes  that  always 

/(a)  =  T*(a)  =  qa  —  r(g)  at  points  a  =  r'(q).  (1-12) 


^Recall  that  we  fixed  a  path  of  Y.  Randomness  is  here  understood  in  choosing  k. 
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Going  through  some  of  the  explicitly  calculated  examples  in  Section  5.5  will  help  de¬ 
mystify  the  Legendre  transform.  A  tutorial  on  the  Legendre  transform  in  contained  in 
Appendix  A  of  [89]. 

From  (1.11)  follows,  that  /(a)  <  f**(a )  =  r*(a)  and  also  that  r(g)  is  a  concave 
function,  hence  continuous  and  almost  everywhere  differentiable. 

1.4  Deterministic  Envelope 

So  far,  all  that  has  been  said  applies  to  any  given  function  or  path  of  a  process.  In  the 
random  case,  one  would  often  like  to  use  a  simple  analytical  approach  in  order  to  gain 
intuition  or  an  estimate  of  /  for  a  typical  path  of  Y. 

To  this  end  we  formulate  a  LDP  for  the  sequence  of  distributions  of  { Iifl> }  where 
randomness  enters  now  through  choosing  k  e  {0, . . .  ,2"  —  1}  as  well  as  through  the 
randomness  of  the  process  itself,  i.e. ,  through  Yt(to)  where  to  lies  in  the  probability 
space  (f2,Pn).  The  moment  generating  function  of  Anfk,to )  =  —  nh^f\u)  ln(2)  with  k 
and  uo  random  is  2_nIEf7[S'(n)(g)].  This  leads  to  defining  the  ‘deterministic  envelope’: 

T(q )  :=liminf — -  log2  IE^S'^g)  (1-13) 

n— xx)  Tl 

and  the  corresponding  ‘rate  function’  F  (see  (3.23)).  As  with  the  pathwise  /(a)  and 
r(g)  we  have  here  again  T(q )  =  F*(q).  More  importantly,  it  is  easy  to  show  that 
r(g,  cu)  >  T(q)  almost  surely  (see  lemma  3.9).  Thus: 

Corollary  1.1.  With  probability  one  the  multifractal  spectra  are  ordered  as  follows:  for 
all  a 

dim(E[a])  <  /(a)  <  r*(a )  <  T*(a),  (1.14) 

provided  that  they  are  all  defined  in  terms  of  the  same  singularity  exponent. 

Great  effort  has  been  spent  on  investigating  under  which  assumptions  equality  holds 
between  some  of  the  spectra,  as  a  matter  of  fact  mostly  between  spectra  based  on  differ¬ 
ent  scaling  exponents.  Indeed,  the  most  interesting  combinations  seem  to  be  dim(P^) 
with  scaling  exponents  h(f'}  and  cqn') ,  and  r*(a)  with  scaling  exponents  vj(fl>  and 
the  former  for  its  importance  in  the  analysis  of  regularity,  the  latter  for  its  numerical 
relevance.  It  has  become  the  accepted  term  in  the  literature  to  say  that  the  multifractal 
formalism  holds  if  any  such  spectra  are  equal;  indeed  they  are  in  a  generic  sense  [52], 
However,  this  terminology  might  sometimes  be  confusing  if  the  nature  of  the  parts  of 
such  an  equality  is  not  indicated.  We  prefer  here  to  call  (1.14)  the  multifractal  formal¬ 
ism:  this  formula  ‘holds’  for  any  fixed  choice  of  a  singularity  exponent  as  is  shown  in 
the  paper. 

1.5  Self-similarity  and  LRD 

The  statistical  self-similarity  as  expressed  in  (1.1)  makes  ffirn,  or  rather  its  increment 
process  a  paradigm  of  long  range  dependence  (LRD).  To  be  more  explicit  let  5  denote 
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some  fixed  lag  and  define  fractional  Gaussian  noise  (fGn)  as 

G(k):=BH((k  +  l)S)~BH(kS).  (1.15) 

Having  the  LRD  property  means  that  the  auto-correlation  ra(k)  :=  IEo[G(n  +  k)G(n)\ 
decays  so  slowly  that  J2krG(k)  =  oo.  The  presence  of  such  strong  dependence  bears 
an  important  consequence  on  the  aggregated  processes 

(fc+l)m— 1 

G{m\k)  :=  -  V  G(i).  (1.16) 

i=km 

They  have  a  much  higher  variance,  and  variability,  than  would  be  the  case  for  a 
short  range  dependent  process.  Indeed,  if  X(k)  are  i.i.d.,  then  X^m\k)  has  variance 
(l/m2)var(A"0+. .  .+ATm_i)  =  (l/m)var(AT).  For  G  we  find,  due  to  (1.1)  and  BH( 0)  =  0, 

=  m2^"2var(G(0)) .  (1.17) 

For  H  >  1/2  this  expression  decays  indeed  much  slower  than  1/m.  As  is  shown  in  [19] 
var(X(m))  ~  m2H  2  is  equivalent  to  rx(k )  ~  k2H  2  and  so,  G{k )  is  indeed  LRD  for 
H  >  1/2  (this  follows  also  directly  from  (7.3)). 

Let  us  demonstrate  with  fGn  how  to  relate  LRD  with  multifractal  analysis  using  only 
that  it  is  a  zero-mean  processes,  not  (1.1).  To  this  end  let  5  =  2-n  denote  the  finest 
resolution  we  will  consider,  and  let  1  be  the  largest.  For  m  =  2*  (0  <  i  <  n )  the  process 
mG^m\k)  becomes  simply  BH((k+  l)m5)  —  BH(km5)  =  BH((k+  l)2*_n)  —  BH{k2l~n). 
But  the  second  moment  of  this  expression  — which  is  also  the  variance —  is  exactly 
what  determines  TQ(2)  (compare  (1.10)).  More  precisely,  using  stationarity  of  G  and 
substituting  m  =  2l,  we  get 

2n-i_i 

2-(n-i)r„( 2)  _  En  [5«-i(2)]  =  ^2  En  [|mG(m)(A;)|2]  =  2n“i22Var  (g (2i))  .  (1.18) 

k= 0 

This  should  be  compared  with  the  definition  of  the  LRD-parameter  H  via 

var(G(m))  ~  m2H~2  or  var(G(2i))  =  2i(2^"2).  (1.19) 

At  this  point  a  conceptual  difficulty  arises.  Multifractal  analysis  is  formulated  in  the 
limit  of  small  scales  (i  — >  — oo)  while  LRD  is  a  property  at  large  scales  (i  — >  oo).  Thus, 
the  two  exponents  H  and  Ta( 2)  can  in  theory  only  be  related  when  assuming  that 
the  scaling  they  represent  is  actually  exact  at  all  scales,  and  not  only  asymptotically. 
When  this  assumption  is  violated,  the  two  approaches  may  provide  strikingly  different 
answers  (compare  Example  7.2). 

In  any  real  world  application,  however,  one  will  determine  both,  H  and  Ta( 2),  by 
Ending  a  scaling  region  i  <  i  <  i  in  which  (1.18)  and  (1.19)  hold  up  to  satisfactory 
precision.  Comparing  the  two  scaling  laws  in  i  yields  Ta( 2)  +  1  —  2  =  2 H  —  2,  or 

Ta( 2)  +  1 


var(G'-m') (()))  =  var  ( — BnimS)]  =  var  ( - Bh(S) 

\m  J  \  m 


2 


(1.20) 
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This  formula  expresses  most  pointedly,  how  multifractal  analysis  goes  beyond  second 
order  statistics :  in  (1.26)  we  compute  with  T(q)  the  scaling  of  all  moments.  The  formula 
(1.20)  is  derived  here  for  zero-mean  processes,  but  can  be  put  on  more  solid  grounds 
using  wavelet  estimators  of  the  LRD  parameter  [4]  which  are  more  robust  than  the  ones 
obtained  through  variance  of  the  increment  process.  The  same  formula  (1.20)  reappears 
also  for  multifractals,  suggesting  that  it  has  some  ‘universal  truth’  to  it,  at  least  in  the 
presence  of  ‘perfect  scaling’  (see  (1.29)  and  (7.25),  but  also  Example  7.2). 

1.6  Multifractal  Processes 

The  most  prominent  examples  where  one  finds  coinciding,  strictly  concave  multifractal 
spectra  are  the  distribution  functions  of  cascade  measures  [64,  56,  15,  33,  6,  82,  49,  91, 
95,  86]  for  which  dim(£dal)  and  T*(a )  are  equal  and  have  the  form  of  a  fl  (see  Figure  6 
and  also  3  (e)).  These  cascades  are  constructed  through  some  multiplicative  iteration 
scheme  such  as  the  binomial  cascade,  which  is  presented  in  detail  in  the  paper  with 
special  emphasis  on  its  wavelet  decomposition.  Having  positive  increments,  however, 
this  class  of  processes  is  sometimes  too  restrictive.  fBm,  as  noted,  has  the  disadvantage 
of  a  poor  multifractal  structure  and  does  not  contribute  to  a  larger  pool  of  stochastic 
processes  with  multifractal  characteristics. 

It  is  also  notable  that  the  first  ‘natural’,  truly  multifractal  stochastic  process  to  be 
identified  was  Levy  motion  [54],  This  example  is  particularly  appealing  since  scaling 
is  not  injected  into  the  model  by  an  iterative  construction  (this  is  what  we  mean  by 
the  term  natural).  However,  its  spectrum  is,  though  it  shows  a  non-trivial  range  of 
singularity  exponents  H(t),  degenerated  in  the  sense  that  it  is  linear. 

Construction  and  Simulation 

With  the  formalism  presented  here,  the  stage  is  set  for  constructing  and  studying  new 
classes  of  truly  multi-fractional  processes.  The  idea,  to  speak  in  Mandelbrot’s  own 
words,  is  inevitable  after  the  fact.  The  ingredients  are  simple:  a  multifractal  ‘time 
warp’,  i.e.,  an  increasing  function  or  process  A4(t)  for  which  the  multifractal  formalism 
is  known  to  hold,  and  a  function  or  process  V  with  strong  mono-fractal  scaling  prop¬ 
erties  such  as  fractional  Brownian  motion  (fBm),  a  Weierstrass  process  or  self-similar 
martingales  such  as  Levy  motion.  One  then  forms  the  compound  process 

V(t)  :=  V(M{t)).  (1.21) 

To  build  an  intuition  let  us  recall  the  method  of  midpoint  displacement  which  can 
be  used  to  define  simple  Brownian  motion  B\i 2  which  we  will  also  call  Wiener  mo¬ 
tion  (WM)  for  a  clear  distinction  from  fBm.  This  method  constructs  Bi/2  iteratively 
at  dyadic  points.  Having  constructed  Bi/2{k2~n)  and  B1/2((k  +  l)2-n)  one  defines 
Bi/2({2k  +  l)2~n_1)  as  (B1/2(k2~n)  +  Bi/2({k  +  l)2_n))/2  +  XkiU.  The  off-sets  Xk)U  are 
independent  zero  mean  Gaussian  variables  with  variance  such  as  to  satisfy  (1.1)  with 
H  =  1/2.  Thus  the  name  of  the  method.  One  way  to  obtain  Wiener  motion  in  multi¬ 
fractal  time  WM(MF)  is  then  to  keep  the  off-set  variables  Xk,n  as  they  are  but  to  apply 


1  Introduction  and  Summary  11 


them  at  the  time  instances  tk,n  defined  by  tk,n  =  M-  1(k2  n),  i.e.,  M.{tk,n )  =  k2  n: 

n  U  \  ^1/2  (4,  n)  +  #1/2  (4+1, n)  ,  ^  ft  nn\ 

&l/2\t2k+l,n+l)  ■—  - ^ - r  Afcjfl. 

This  amounts  to  a  randomly  located  random  displacement ,  the  location  being  deter¬ 
mined  by  Ad.  Indeed,  (1.21)  is  nothing  but  a  time  warp. 

An  alternative  construction  of  ‘warped  Wiener  motion’  WM(MF)  which  yields  an 
equally  spaced  sampling  as  opposed  to  the  samples  #1/2  (4, n)  provided  by  (1.22)  is 
desirable.  To  this  end,  note  first  that  the  increments  of  WM(MF)  become  independent 
Gaussians  once  the  path  of  A4(t)  is  realized.  To  be  more  precise,  fix  n  and  let 

Q{k)  :=  B((k  +  l)2-n)  -  B{k2-n )  =  B1/2(M(k  +  l)2~n))  -  B1/2(M(k2~n)).  (1.23) 

For  a  sample  path  of  Q  one  starts  by  producing  first  the  random  variables  M.(k2~n). 
Once  this  is  done,  the  Q{k)  simply  are  independent  zero- mean  Gaussian  variables  with 
variance  \M.(k  +  l)2-n))  —  M.{k2~n)\.  This  procedure  has  been  used  in  Figure  3. 


Global  analysis 

For  the  right  hand  side  (RHS)  of  the  multifractal  formalism  (1.14),  we  need  only  to 
know  that  V  is  an  H- sssi  process,  meaning  that  the  increment  V(t  +  u)  —  V (t)  is  equal 
in  distribution  to  uHV(  1)  (compare  (1.1)).  Assuming  independence  between  V  and  A4 
a  simple  calculation  reads  as 


2n  — 1 


En  \V((k  +  I)2”")  -  V(k2 

k= 0 
2n  — 1 

EE  \\V(M({k  +  l)2-n))-V(M{k2~n))\q  M(k2~n),  M((k  +  l)2"n) 

k= 0 
2n—l 

E  [I M((k  +  l)2_n)  -  M{k2~n)\qH]  IE  [^(l)!9] .  (1.24) 


k= 0 


Here,  we  dealt  with  increments  |V((fc  +  l)2  n)  —  V{k2  n)|  for  the  ease  of  notation.  With 

,  (n)  (n) 

little  more  effort  they  can  be  replaced  by  suprema,  i.e.,  by  2~nnk  s  or  even  by  2~nwk 
for  certain  wavelet  coefficients  and  under  appropriate  assumptions  (see  theorem  8.5). 
It  follows,  e.g.,  for  h^\  that 


Warped  H- sssi:  Thy(q) 


ThM^n)  if  En  [i  supo^!  V(t)\q]  <  oo 
— oo  else. 


(1.25) 


Simple  H- sssi  process:  When  choosing  the  deterministic  warp  time  M(t)  =  t  we 
have  TM(q)  =  q  —  1  since  S^(q)  =  const2”  •  2~nq  for  all  n.  Also,  V  —  V.  We  obtain 
T,vf  (qH)  =  qH  —  1  which  has  to  be  inserted  into  (1.25)  to  obtain 


(  qH  -1  if  Eq  [j  sup,,^!  V (t) \q] 
\  —oo  else. 


Simple  H- sssi: 


<  (X) 


(1.26) 
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Local  analysis  of  warped  fBm 

Let  us  now  turn  to  the  special  case  where  V  is  fBm.  Then,  we  use  the  term  FB(MF)  to 
abbreviate  fractional  Brownian  motion  in  multifractal  time:  B(t)  =  BH(A4(t)).  First, 
to  obtain  an  intuition  on  what  to  expect  from  the  spectra  of  B  let  us  note  that  the 
moments  appearing  in  (1.25)  are  finite  for  all  q  as  we  will  see  in  lemma  7.4.  Applying 
the  Legendre  transform  yields  easily  that 

W  =  inf (qa  -  TB(q ))  =  inf(go  -  TM{qH))  =  T*M{a/H\  (1.27) 

q  q 

which  is  valid  for  all  a  €  IR  for  which  the  second  equality  holds,  i.e.,  for  which  the 
inbmum  is  attained  for  q  values  in  the  range  where  TB(q )  is  finite.  In  particular,  for 
Brownian  motion  (fBm  with  H  =  1/2)  it  holds  for  all  a  (compare  lemma  7.4). 


(a) 


(b) 


(d) 


time  lag 


a 


FIGURE  3.  Left:  Simulation  of  Brownian  motion  in  binomial  time  (a)  Sampling  of 
Mb{(k  +  l)2_n)  —  Mb{k2~n)  (k  =  0,...,2n  —  1),  indicating  distortion  of  dyadic 

time  intervals  (b)  Mb((k2~n)):  the  time  warp  (c)  Brownian  motion  warped  with  (b): 
B(k  2~n)  =  B1/2(Mb(k2~n )) 

Right:  Estimation  of  dim(Eg^)  via  r*,  B  (d)  Empirical  correlation  of  the  Haar  wavelet  coeffi¬ 
cients.  (e)  Dot-dashed:  (from  theory),  dashed:  Tg(a)  =  (a/ H)  Solid:  the  estimator 

Tw  B  obtained  from  (c).  (Reproduced  from  [40].) 

Second,  towards  the  local  analysis  we  recall  the  uniform  and  strict  Holder  continuity 
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of  the  paths  of  fBm.  In  theorem  7.3  we  state  a  precise  result  due  to  Adler  [5]  which 
reads  roughly  as 

sup  | B(t  +  u)  —  B(t) |  =  sup  | +  u ))  —  BH(M(t))\  ~  sup  \M{t  +  u)  —  M{t)\H . 

|u|<£  |u|<5  1^1^^ 

This  is  the  key  to  conclude  that  BH  simply  squeezes  the  Holder  regularity  exponents 
by  a  factor  H .  Thus, 

hB{t)  =  H-hM{t),  E[tiH]  =  E[*\ 

and,  consequently,  analogous  to  (1.27), 


dim(ifj^)  =  dim 


Figure  3  (d)-(e)  displays  an  estimation  of  dim(4^)  using  wavelets  which  agrees  very 
closely  with  the  form  dim (E^H^)  predicted  by  theory.  (For  statistics  on  this  estimator 
see  [40,  41].) 

Combining  this  with  corollary  1.1  and  (1.27)  we  may  conclude: 

Corollary  1.2  (Fractional  Brownian  Motion  in  Multifractal  Time). 

Let  Bh  denote  fBm  of  Hurst  parameter  H.  Let  M.{t)  be  of  almost  surely  continuous 
paths  and  independent  of  BH.  Set  Bit )  =  BH(A4(t))  and  consider  a  multifractal  anal¬ 
ysis  using  tiff'1 .  Then,  the  multifractal  warp  formalism 


dim(4a])  =  fB(a)  =  Tft{a)  =  T*(a)  =  T*M{a/H) 


(1.28) 


holds  for  any  path  and  any  a  for  which  dim (44)  =  Tf^ia/  H)  =  Tg(a). 

The  condition  on  a  ensures  that  equality  holds  in  the  multifractal  formalism  for  Jvl 
and  that  the  relevant  moments  are  finite  so  that  (1.27)  holds.  If  satisfied,  then  the 
local,  or  fine,  multifractal  structure  of  B  captured  in  dim(£’^)  on  the  left  side  in  (1.28) 
can  be  estimated  through  grain  based,  simpler  and  numerically  more  robust  spectra 
on  the  right  side,  such  as  r£(a)  (compare  Figure  3  (e)). 

Moreover,  the  ‘warp  formula’  (1.28)  is  appealing  since  it  allows  to  separate  the  LRD 
parameter  of  fBm  and  the  multifractal  spectrum  of  the  time  change  A4.  Indeed,  pro¬ 
vided  that  M.  is  almost  surely  increasing  one  has  Tm(1)  =  0  since  S'^(O)  =  -M(l) 
for  all  n.  Thus,  T^(l/H)  =  0  exposes  the  value  of  H.  Alternatively,  the  tangent  at 
through  the  origin  has  slope  1  / H.  Once  H  is  known  follows  easily  from  T$. 

Simple  fBm:  When  choosing  the  deterministic  warp  time  A4 (t)  =  t  we  have  B  =  Bh 
and  TBH(q)  =  qH  —  1  as  in  (1.26).  In  the  special  case  of  Brownian  motion  (H  =  1/2)  we 
may  apply  (1.28)  for  all  a  showing  that  all  /q/Mjased  spectra  consist  of  the  point  (H,  1) 
only.  This  makes  the  mono- fractal  character  of  this  process  most  explicit.  In  general, 
however,  artifacts  which  are  due  mainly  to  diverging  moments  may  distort  this  simple 
picture  (see  Section  7.3). 
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LRD  and  estimation  of  warped  Brownian  motion 


Let  Q[k)  :=  B([k+l)2~n)—B{k2~n)  be  ‘fGn  in  multifractal  time’  (see  (1.23)  for  the  case 
H  =  1/2).  Calculating  auto-correlations  explicitly,  lemma  8.8  shows  that  G  is  second 
order  stationary  provided  At  has  stationary  increments.  Assuming  IE[.Ad(<s)2^]  =  const- 
gT(2H)+ the  correlation  of  Q  is  of  the  form  of  ordinary  fGn,  but  decaying  as  rg(k)  ~ 
k2Hg~ 2  where 


Hg  = 


Tm(2H)  +  1 
2 


(1.29) 


Let  us  discuss  some  special  cases.  An  obvious  choice  for  a  subordinator  At  is  Levy  mo¬ 
tion,  an  //'-self-similar,  1/LF-stable  process.  It  has  independent,  stationary  increments. 
Since  the  relation  (1.1)  holds  with  H'  as  the  scaling  parameter,  we  have  T(q)  =  qH'  —  1 
from  (1.26).  Moreover,  Ai(s)2H  is  equal  in  distribution  to  (sH  Ai{l))2H  and  indeed 
!E[A/f(s)2fl]  =  const  •  s2HH'  =  const  •  sTt2H)+1  This  expression  is  finite  for  2 H  <  1  / H' . 
In  summary,  Hg  =  II IV  <1/2. 

For  a  continuous,  increasing  warp  time  At,  on  the  other  hand,  we  have  always 
Tvf  (0)  =  —1  and  TM(\)  =  0.  (Levy  motion  is  discontinuous;  it  is  increasing  for  H'  <  1, 
in  which  case  T(l)  is  not  defined.)  Exploiting  the  concave  shape  of  we  find  that 
H  <  Hg  <  1/2  for  0  <  H  <  1/2,  and  1/2  <  Hg  <  H  for  the  LRD  case  1/2  <  H  <  1. 

Especially  in  the  case  of  H  =  1/2  (‘white  noise  in  multifractal  time’)  G(k)  becomes 
uncorrelated  (see  also  (8.20)).  Notably,  this  is  a  stronger  statement  than  the  observation 
that  the  Gik)  are  independent  conditioned  on  At  (compare  Section  1.6).  As  a  particular 
consequence,  wavelet  coefficients  will  decorrelate  fast  for  the  compound  process  G,  not 
only  when  conditioning  on  At  (compare  Figure  3  (d)).  This  is  favorable  for  estimation 
purposes  as  it  reduces  the  error  variance.  Finally,  for  increasing  A4  we  have  T(l)  =  0 
and  the  requirements  for  (1.29)  reduce  to  the  simple  IE[Al(s)]  =  s.  Multiplicative 
processes  with  this  property  (as  well  as  stationary  increments)  have  been  introduced 
recently  [14,  70,  74,  105]. 

Though  seemingly  obvious  it  should  be  pointed  out  that  the  vanishing  correlations 
of  G  in  the  case  H  =  1/2  should  not  be  taken  as  an  indication  of  independence.  After 
all,  G  becomes  Gaussian  only  when  conditioning  on  knowing  At.  A  strong,  higher  order 
dependence  in  G  is  hidden  in  the  dependence  of  the  increments  of  At  which  determine 
the  variance  of  G(k)  as  in  (1.23).  Indeed,  Figure  3  (c)  shows  clear  phases  of  monotony 
of  B  indicating  positive  dependence  in  its  increments  G,  despite  vanishing  correlations. 
Mandelbrot  calls  this  the  ‘blind  spot  of  spectral  analysis’. 


Multifractals  in  multifractal  time 

Despite  of  its  simplicity  the  presented  scheme  for  constructing  multifractal  processes 
allows  for  various  play-forms  some  of  which  arc  little  explored.  First  of  all,  for  simulation 
purposes  one  might  subject  the  randomized  Weierstrass- Mandelbrot  function  to  time 
change  rather  than  ffim  itself. 

Next,  we  may  choose  to  replace  ffim  with  a  more  general  self-similar  process  (7.1) 
such  as  Levy  motion.  Difficulties  arise  here  since  Levy  motion  is  itself  a  multifractal 
with  varying  Holder  regularity  and  the  challenge  lies  in  studying  which  exponents  of 
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the  ‘multifractal  time’  and  the  motion  are  most  likely  to  meet.  A  solution  for  the 
spectrum  /(a)  is  given  in  corollary  8.13  which  actually  applies  to  arbitrary  processes  Y 
with  stationary  increments  (compare  theorem  8.15)  replacing  fBm.  In  its  most  compact 
form  onr  result  reads  as: 

Corollary  1.3  (Levy  motion  in  multifractal  time).  Let  Lh  denote  Levy  stable 
motion  and  let  A4  be  a  binomial  cascade  (see  5.1)  independent  of  LH  and  set  V(t)  = 
Lu^M-it)).  Then,  for  almost  all  paths 

fv(a)  =  Tv(a)  =  TZ(a )  (1.30) 


for  all  a  where  Ty  >  0.  The  envelope  Ty  can  be  computed  through  the  warp  formula 


Tv(q)  =  Tm  (TLH{q)  +  l)  • 


(1.31) 


Recall  (1.26)  for  a  formula  of  TjJfI ,  which  is  generalized  in  (7.10).  As  the  paper  shows 
(1.30)  and  (1.31)  hold  actually  in  more  generality. 

Finally,  for  y(t)  =  Y(A4(t))  where  Y  and  M.  are  both  almost  surely  increasing,  i.e., 
multifractals  in  the  classical  sense,  a  close  connection  to  the  so-called  ‘relative  multi¬ 
fractal  analysis’  [95]  can  be  established  using  the  concept  of  inverse  multifractals  [94]: 
understanding  the  multifractal  structure  of  y  is  equivalent  to  knowing  the  multifractal 
spectra  of  Y  with  respect  to  AW1,  the  inverse  function  of  M..  We  will  show  how  this 
can  be  resolved  in  the  simple  case  of  binomial  cascades. 


2  Singularity  Exponents 

For  simplicity  we  consider  processes  Y  over  a  probability  space  (0,  T ,  Pq)  and  defined 
on  a  compact  interval,  which  we  assume  without  loss  of  generality  to  be  [0,1].  Gen¬ 
eralization  to  higher  dimensions  is  straightforward  and  extending  to  processes  defined 
on  IR  is  simple  and  will  be  indicated. 

2.1  Holder  Continuity 

As  pointed  out  in  the  introduction,  the  erratic  behavior  of  a  continuous  process  Y(t) 
maybe  indicative  of  crucial  properties  with  relevance  in  various  applications.  This  local 
behavior  of  Y  at  a  given  time  t  can  be  characterized  to  a  first  approximation  by 
comparison  with  an  algebraic  function  as  follows: 

Definition  2.1.  A  function  or  the  path  of  a  process  Y  is  said  to  be  in  C1)  if  there  is  a 
polynomial  Pt  such  that 

\Y(u)  —  Pt{u)\  <  C\u  —  t\h 

for  u  sufficiently  close  to  t.  Then,  the  degree  of  local  Holder  regularity  of  Y  at  t  is 

H(t)  :=  sup {h  :  Y  G  Cf}.  (2.1) 
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As  usual,  let  |_Al  denote  the  largest  integer  smaller  or  equal  to  x.  If  the  Taylor 
polynomial  of  degree  [H{t)\  exists,  then  P  is  necessarily  that  Taylor  polynomial.  As 
the  example  Y(t)  =  1  + 1  + 12  +  f3'5  sin(l/f)  shows,  P  may  be  different  from  the  Taylor 
polynomial  if  Y  does  not  have  sufficient  degree  of  smoothness.  Here,  Y(t)  has  only 
one  derivative  at  t  —  0,  and  its  Taylor  polynomial  at  t  =  0  is  u  i— >  1  +  u  while 
Po(u)  =  1  +  u  +  u2. 

Of  special  interest  for  our  purpose  is  the  case  when  the  approximating  polynomial 
Pt  is  a  constant,  i.e.,  Pt{u)  =  Y(t),  in  which  case  H(t )  can  be  computed  easily.  To  this 
end: 

Definition  2.2.  Let  us  agree  on  the  convention  log(0)  =  —  oo  and  set 

h(t)  liminf  - - log2  sup  \Y(u)-Y(t)\.  (2.2) 

log2(2e)  |u-t|<e 

Note  first  that  for  any  h  <  h(t )  we  have  | Y(u)  —  Y(t) \  <  C\u  —  t\h,  thus  Y  e  C f. 
We  conclude  that  always 

h(t)  <  H(t).  (2.3) 

It  is  easy  to  see  that  if  Pt  is  known  to  be  a  constant  then  H(t)  =  h(t).  On  the  other 
hand,  we  show: 

Lemma  2.3.  If  h(t)  £  IN  then  Pt  is  a  constant,  and  hit )  =  Hit). 

As  the  example  Y(t)  =  t 2  +  t2A  shows  the  conclusion  does  not  necessarily  hold  when 
h(t)  e  IN. 

Proof 

We  will  show  the  dual  statement:  If  H{t)  >  h(t)  then  h{t)  must  be  an  integer.  Together 
with  (2.3)  this  will  certainly  establish  the  lemma. 

So,  we  assume  H(t)  >  h(t).  Then  there  is  h  >  h(t )  and  a  polynomial  Pt(-)  such  that 
| Y(u)  —  Ptiu)  |  <  C\u  —  t\h.  Note  that  Pt  cannot  be  a  constant:  if  it  were  constant,  then 
h  <  h(t)  due  to  the  very  definition  of  h(t).  Thus,  we  may  write  the  error-minimizing 
polynomial  as  Ptfu)  =  Y(t)+(u— t)m-Q(u)  for  some  integer  m  >  1  and  some  polynomial 
Q  without  zero  at  t.  Assume  first  that  m  <  lift)  and  choose  h'  such  that  m  <  h!  <  hit). 
Writing  Y(u)  —  Ptfu)  =  (' Y(u )  —  Y(t ))  —  (Ptfu)  —  Y{t)),  the  first  term  is  smaller  than 
| u  —  t\h  and  the  second  term,  decaying  as  C\u  —  t\m,  governs.  Whence  h  —  rn  <  h(t), 
a  contradiction  against  the  assumption  h  >  h{t).  Assuming  m  >  hit)  choose  h'  such 
that  m  >  h!  >  hit)  and  a  sequence  un  such  that  \Y(un)  —  Y(t)  \  >  \un  —  t\h  .  Then, 

| Y(un)  —  Ptfun) |  >  (1/2)| un  —  t\h'  for  large  n  and  h  <  h! .  Letting  h!  — >  hit)  we  get 
again  a  contradiction.  We  conclude  that  hit)  equals  m.  <0 

An  essential  simplification  for  both,  analytical  and  empirical  study,  is  to  replace  the 
continuous  limit  in  (2.2)  by  a  discrete  one.  To  this  end  we  introduce  some  notation 

Definition  2.4.  Let  kn(t)  :=  [t2n\.  Then,  kn(t)  is  the  unique  integer  such  that 

t  e  /'”>  :=  \kn2-",  (kn  +  l)2-"[. 


(2.4) 


2  Singularity  Exponents  17 


As  n  increases  the  intervals  lj"f  form  a  nested  decreasing  sequence  (compare  Fig¬ 
ure  4).  Now,  when  defining  a  discrete  approximation  to  hit)  we  have  to  imitate  in 
a  discrete  manner  a  ball  around  t  over  which  we  will  consider  the  increments  of  Y. 
Accounting  for  the  fact  that  t  could  lay  very  close  to  the  boundary  of  ij^,  we  set 

Definition  2.5.  The  coarse  Holder  exponents  ofY  are 

h^  :=  -ilog2sup{|y(u)  -  Y(t)\  :  u  €  [(kn  -  1)2“",  (kn  +  2)2“"[}.  (2.5) 

To  compare  the  limiting  behavior  of  these  exponents  with  h{t)  we  choose  n  such  that 
2~n+1  <  e  <  2~n+2.  We  have 


[{kn  -  1)2“",  (kn  +  2)2“"[c  [t  +  e,t-  e[c  [(fc„_ 2  -  1)2“"+2,  ( kn_2  +  2)2“n+2[. 

from  which  it  follows  immediately  that 

Lemma  2.6. 

h(t)  =  lim  inf  h ^ 

n—>o o  n 

It  is  essential  to  note  that  the  countable  set  of  numbers  lif]  contains  all  the  scaling 
information  of  interest  to  us.  Being  defined  pathwise,  they  are  random  variables. 

2.2  Scaling  of  Wavelet  Coefficients 

The  discrete  wavelet  transform  represents  a  1-d  process  Y(t)  in  terms  of  shifted  and 
dilated  versions  of  a  prototype  bandpass  wavelet  function  -0(f),  and  shifted  versions  of 
a  low-pass  scaling  function  0(f)  [23,  106].  Made  precise  in  the  vocabulary  of  Hilbert 
spaces:  For  special  choices  of  the  wavelet  and  scaling  functions,  the  atoms 

tMt)  :=  2j'2if{2jt  -  k)  ,  fc(t)  :=  2j'2  4>{2H  -  k)  ,  j,ke  Z  (2.6) 


form  an  orthonormal  basis  and  we  have  the 

Y  (t)  =  w 

representations  [23,  106] 

OO 

+  'C'j,klfj,k{t), 

(2.7) 

k 

II 

with 

r  n 

Cjlk  :=  J  Y(t)  if  j, kit)  df, 

Dj,k  :=  J  Y(t)  4>j,k(t)  df. 

(2,8) 

For  a  wavelet  0(f)  centered  at  time  zero  and  frequency  /o,  the  wavelet  coefficient 
Cj^k  measures  the  signal  content  around  time  2 “J/c  and  frequency  2J,/0.  The  scaling 
coefficient  D]  k  measures  the  local  mean  around  time  2 “J/c.  In  the  wavelet  transform, 
j  indexes  the  scale  of  analysis:  Jq  can  be  chosen  freely  and  indicates  the  coarsest  scale 
or  lowest  resolution  available  in  the  representation. 

Compactly  supported  wavelets  are  of  especial  interest  in  this  paper  (see  [23].  The 
Haar  scaling  and  wavelet  functions  (see  Figure  4(a))  provide  the  simplest  example  of 
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such  an  orthonormal  wavelet  basis:  0  is  the  indicator  function  of  the  unit  interval,  while 
0(f)  =  0(2t)  —  0(2t  —  1).  For  a  process  supported  on  the  unit  interval  one  may,  thus, 
choose  Jq  =  0.  The  supports  of  the  fine-scale  scaling  functions  nest  inside  the  supports 
of  those  at  coarser  scales;  this  can  be  neatly  represented  by  the  binary  tree  structure  of 
Figure  4(b).  Row  (scale)  j  of  this  scaling  coefficient  tree  contains  an  approximation  to 
Y{t)  of  resolution  2~f  Row  j  of  the  complementary  wavelet  coefficient  tree  (not  shown) 
contains  the  details  in  scale  j  +  1  of  the  scaling  coefficient  tree  that  are  suppressed  in 
scale  j.  In  fact,  for  the  Haar  wavelet  we  have 


Dj,k  —  2  1l2(yDj+i£k  +  Dj+i,2k+i), 
Cj,k  =  ^  l^2{Dj+it2k  ~  Dj+ l,2fc+l)- 


(2.9) 


2 


j/2 


0 


k2'J  (k+l)2 


2/2 

... 

n 

0 

k2'j 

□ 

(k+l)2j 

(a) 


FIGURE  4.  (a)  The  Haar  scaling  and  wavelet  functions  <fj,k(t)  and  'if3y{t).  (b)  Binary  tree 
of  scaling  coefficients  from  coarse  to  fine  scales. 


Wavelet  decompositions  contain  considerable  information  on  the  singularities  of  a 
process  Y.  Indeed,  adapting  the  argument  of  [51,  p.  291]  (note  the  L2  wavelet  normal¬ 
ization  used  here  —  as  opposed  to  L 1  in  [51]  )  we  find 

Lemma  2.7.  Fix  t  and,  let  kn  =  kn(t )  as  in  (2.4).  If  l'K(s)  —  Y(t) \  =  0(|s  —  t\h)  for 
s  — >  t,  and  if  if  is  a  compactly  supported  function  with  J  if  =  0  and  f  \if\  <  oo,  then 


2n/2 1  a 


n,kn 


=  2r 


Y(s)  if(2ns  —  k)  ds 


=  0( 2~nh)  for  n  oo. 


(2.10) 


Proof 

The  compact  support  of  if  is  assumed  only  for  simplicity  of  the  argument  (compare 
[51]),  and  we  take  it  to  be  [0, 1].  Then,  if{2ns  —  k)  =  0  for  s  <£  ljfl) .  Also,  \s  — 1|  <  2~n 
for  all  s  G  l['‘\  These  facts,  together  with  f  Y (t)if(s)ds  =  0  and  the  given  estimate  on 
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Y  allow  to  conclude  as 


2”n/2 1 C, 


n,kn 


An) 


(Y(s)  —  Y(t))  if(2ns  —  k)  ds 


<  C- 


=  C  ■  2 


An) 

Lk 

—nh 


|s  —  t\h  \f>(2ns  —  k)\ds 


T(n ) 


\if{2ns  -  k)  |  ds  =  C  •  2 


—nhrx—n 


(s)  |  ds 


0 

Following  the  proof  of  [51,  p.  291]  more  closely  it  is  easily  seen  that  the  assumption 
of  xf  being  compactly  supported  is  not  really  needed.  Indeed,  a  fast  decay  is  sufficient 
and  the  result  holds  also  for  functions  which  don’t  necessarily  form  a  basis  such  as 
derivatives  of  the  Gaussian  exp(— x2).  To  distinguish  such  functions  from  the  orthogonal 
wavelets  we  will  address  them  as  ‘analyzing  wavelets’. 

In  order  to  invert  lemma  2.7  and  infer  the  Holder  regularity  of  Y  from  the  decay  of 
wavelet  coefficients  one  needs  the  representation  (2.7),  sufficient  wavelet  regularity  as 
well  as  some  knowledge  on  the  decay  of  the  maximum  of  the  wavelet  coefficients  in  the 
vicinity  of  t,  as  developed  in  the  pioneering  work  of  [78]  (see  also  [51]  and  [23,  Thm. 
9.2]). 

All  this  suggests  that  the  left  hand  side  of  (2.10)  could  produce  an  alternative  useful 
description  of  the  local  behavior  oscillatory  behavior  of  Y. 

Definition  2.8.  The  coarse  wavelet  singularity  exponents  of  Y  are 

WK.  ~“l°g2  |2n/2C'ri)fen|  •  (2.11) 


The  local  singularity  exponent  of  wavelet  coefficients  is  then  w(t ) 


lim  inf  w 

77^00 


(n) 

kn 


Indeed,  while  the  coarse  Holder  exponents  give  the  exact  Holder  regularity  but 
only  under  the  assumption  of  constant  approximating  polynomials,  using  wavelets  has 
the  advantage  of  yielding  an  analysis  which  is  largely  unaffected  by  polynomial  trends 
in  Y.  This  is  due  to  vanishing  moments  f  tm'ip{t)dt  =  0  which  are  typically  built  into 
wavelets  [23].  However,  for  a  reliable  estimation  of  true  Holder  continuity  through 
wavelets  one  has  to  employ  the  lines  of  maxima,  a  method  pioneered  convincingly  in 
[78]  (see  also  [51,  23,  53,  7]). 

In  any  case,  the  decay  of  wavelet  coefficients  is  interesting  in  itself  as  it  relates  to 
LR.D  (compare  [4]  and  Section  7.4)  and  regularity  spaces  such  as  Besov  spaces  [89,  30]. 


2.3  Other  Singularity  Exponents 

The  ‘classical’  multifractal  analysis  of  a  singular  measure  /i  on  [0,1],  translated  into 
the  notations  used  here,  has  always  been  concerned  with  the  study  of  the  singularity 
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structure  of  its  primitive 

M{t)=  [  Mds)  =  /i([0,t]),  (2.12) 

Jo 

which  is  an  almost  surely  increasing  process.  Consequently,  the  supremum  in  sim¬ 
plifies  to  | Ad  ((kn  +  2)2~n)  —  A4((kn  —  l)2_n)|  and  one  is  lead  to  introduce  yet  another 
singularity  exponent: 

Definition  2.9.  The  coarse  increment  exponents  of  A4  are 

4?  :=  -^log2|Af  ((kn  +  l)2~n)  -  M(kn2~n)\  =  -4  log. :2/x(/£}).  (2.13) 


The  local  singularity  exponent  of  increments  is  then  aft ) 


lim  inf  a 

n— >oo 


(n) 

kn 


As  we  will  elaborate  in  lemma  5.5  the  increment  and  Holder  exponents  off1  and  h 4 
provide  largely  the  same  analysis  for  increasing  Ad,  however,  there  is  a  crucial  lesson 
to  be  learned  for  measures  /i  with  fractal  support  in  Section  5.6. 

For  further  examples  of  singularity  exponents  we  would  like  to  refer  to  [63]  which 
treats  the  case  of  exponents  which  are  so-called  Choquet  capacities,  a  notion  which  is 
not  needed  to  develop  the  multifractal  formalism  as  we  will  show  in  this  paper. 

Also,  [87]  considers  an  arbitrary  function  £(/)  from  the  space  of  all  intervals  to 
1R+  (instead  of  only  the  ijf1'1 )  and  develops  a  multifractal  formalism  similar  to  the  one 
presented  here.  There,  it  is  suggested  to  consider  the  oscillations  of  Y  around  the  mean, 


i.e., 


fjY(s)ds 

I/I 


df 


(2.14) 


This  gives  raise  to  the  singularity  exponent  —  (1  /n)  log2(£(/4))  which  is  of  particular 
interest  since  it  can  be  used  to  define  oscillation  spaces  such  as  Sobolev  spaces  and 
Besov  spaces.  Alternatively,  interpolating  Y  in  the  interval  /  by  the  linear  function 
a i  +  bit,  one  could  use 


(Yft)  -  ( a i  +  bit))2  df 


1/2 


(2.15) 


This  £(/)  measures  the  variability  of  Y  and  is  related  to  the  dimension  of  the  paths  of 
Y.  In  the  definitions  (2.14)  and  (2.15)  constant,  resp.  linear  terms  are  subtracted  from 
Y.  This  may  remind  one  of  the  use  of  wavelets  with  one,  resp.  two  vanishing  moments. 

In  conclusion,  there  are  various  useful  notions  of  singularity  exponents  which  may 
provide  a  characterization  of  a  process  Y  of  relevance  in  particular  applications.  Being 
aware  that  these  descriptions  may  very  well  differ  for  the  same  Y  according  to  one’s 
choice  of  an  exponent,  we  do  not  attempt  to  value  one  over  the  other,  but  rather  present 
some  aspects  of  multifractal  analysis  which  are  valid  for  any  such  choice,  i.e.,  a  form 
of  the  celebrated  multifractal  formalism. 
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3  Multifractal  Analysis 

Multifractal  analysis  has  been  discovered  in  the  context  of  fully  developed  turbulence 
[64,  36,  42,  46,  45]  and  subsequently  further  developed  in  physical  and  mathematical 
circles  (see  [15,  56,  13,  12,  16,  102,  9,  66,  49,  78,  94,  22,  58,  33,  6,  53,  82,  25,  86,  8,  54,  31, 
74]  to  give  only  a  short  list  of  some  relevant  work  done  in  this  area).  At  the  beginning 
stands  the  discovery  that  on  fractals  local  scaling  behavior  as  measured  by  exponents 
h^\  ay1"1  or  w^\  is  not  uniform  in  general.  In  other  words,  h(t),  aft )  and  w(t)  are 
typically  not  constant  in  t  but  assume  a  whole  range  of  values,  thus  imprinting  a  rich 
structure  on  the  object  of  interest  [69,  80,  2,  88].  This  structure  can  be  characterized 
either  in  geometrical  terms  making  use  of  the  concept  of  dimension,  or  in  statistical 
terms  based  on  sample  moments.  A  tight  connection  between  these  two  descriptions 
emerges  from  the  multifractal  formalism. 

As  we  will  see,  as  far  as  the  validity  of  the  multifractal  formalism  is  concerned  there 
is  no  restriction  in  choosing  a  singularity  exponent  which  seem  fit  for  describing  scaling 
behavior  of  interest,  as  long  as  one  is  consistent  in  using  the  same  exponents  for  both, 
the  geometrical  and  the  statistical  description.  To  express  this  fact  we  consider  in  this 
section  the  arbitrary  coarse  singularity  exponent 

4n)  (fc  =  0,...,2"-l,  nelN),  (3.1) 

which  may  be  any  sequence  of  random  variables.  To  keep  a  connection  with  what  was 
said  before  think  of  ,s,[n'1  as  representing  a  coarse  singularity  exponent  related  to  the 
oscillations  of  Y  over  the  dyadic  interval  ijf1 .  To  accommodate  processes  which  are 
constant  over  some  intervals  we  explicitly  allow  sjy  to  take  the  value  oo. 

3.1  Dimension  based  Spectra 

The  strongest  interests  of  the  mathematical  community  are  in  the  various  measure 
theoretical  dimensions  of  sets  A  A  which  are  defined  pathwise  in  terms  of  limiting 
behavior  of  s[]^  as  n  — >  oo,  as  follows 

A^  \=  {t  :  liminf  sA  =  a}, 

n—>oo  n 

rH  {£  .  ijm  An)  —  aj  (3.2) 

n— >o o  n 

These  sets  are  typically  ‘fractal’  meaning  loosely  that  they  have  a  complicated  geo¬ 
metric  structure  and  more  precisely  that  their  dimensions  are  non-integer.  A  compact 
description  of  the  singularity  structure  of  Y  is,  therefore,  in  terms  of  the  dimensions  of 
A  A  and  A"  A. 

Definition  3.1.  The  Hausdorff  spectrum  is  the  function 

a  ^  dim(AW),  (3.3) 

where  dirn(A)  denotes  the  Hausdorff  dimension  of  the  set  A  [103]. 
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The  sets  (a  G  IR)  form  a  multifractal  decomposition  of  the  support  of  Y,  i.e., 
they  are  disjoint  and  their  union  is  the  support  of  Y.  We  will  loosely  address  Y  as 
a  multifractal  if  this  decomposition  is  rich,  i.e.  if  the  sets  (a  G  IR)  are  highly 
interwoven  or  even  dense  in  the  support  of  Y. 

We  should  point  out  that  the  study  of  singular  measures  (deterministic  and  random) 
has  often  focussed  on  the  simpler  sets  and  their  spectrum  dim(iC^)  [56,  15,  33, 
6,  82,  94,  91,  93,  8].  However,  lemma  3.3  which  is  established  below  allows  to  extend 
most  of  these  results  in  order  to  provide  formulas  for  dim(i^al)  also. 


3.2  Grain  based  Spectra 

An  alternative  to  the  above  geometrical  description  of  the  singularity  structure  relies 
on  counting:* 

N(n\a,  e)  :=  4f{k  —  0, . . . ,  2n  —  1  :  a  —  £  <  <  a  +  e}.  (3.4) 

Note,  that  infinite  have  no  influence  here.  Indeed,  computing  multifractal  spectra 
at  o  =  oo  requires  usually  special  attention  (see  [94]). 

Definition  3.2.  The  grain  based  spectrum  is  the  function 

a)  :=  limlimsup  —  log2  N^n\a,  e).  (3. 

NO  n — >oo  U 

To  establish  some  of  the  almost  sure  pathwise  properties  it  is  convenient  to  introduce 
also 

/(a)  :=  limliminf  —  log2  N^fa,  e)  (3.6) 

—  frj.0  n— >  oo  fl 


This  approach  has  grown  out  of  the  difficulties  involved  with  computation  of  Haus- 
dorff  dimensions,  in  particular  in  any  real  world  applications.  Using  a  mesh  of  given 
grain  size  as  in  (3.4)  instead  of  arbitrary  coverings  as  in  dim(TJ[“l)  leads  generally  to 
more  simple  notions.  However,  /  should  not  be  regarded  as  an  auxiliary  vehicle  but 
rather  meriting  attention  in  its  own  right.  This  point  was  already  made  in  Section  1.2, 
and  we  hope  to  make  it  stronger  as  we  proceed  in  our  presentation. 

Our  first  remark  on  /(a)  concerns  the  fact  that  the  counting  used  in  its  definition, 
i.e.,  N^n\a,  e)  may  be  used  to  estimate  box  dimensions.  Based  on  this  fact  it  was  shown 
in  [94]  that 

dim (A'W)  <f(a).  (3.7) 

Note  that  the  set  in  (3.7)  is  a  subset  of  E^  since  in  its  points  the  sequence  of  sf>}  is 
actually  required  to  converge.  Also,  we  will  later  need  a  lower  bound  on  /(a).  Therefore, 
we  provide  two  formulas  that  are  sharper  than  (3.7). 


fMore  generally,  using  c-ary  intervals  in  Euclidean  space  IRc/  kn  will  range  from  0  to  cnd  —  1. 
Logarithms  will  have  to  be  taken  to  the  base  c  since  we  seek  the  asymptotics  of  j\Rn)  (a,  e)  in 
terms  of  a  powerlaw  of  resolution  at  stage  n,  i.e.,  N^n\a:e)  ~  cnRab  The  maximum  value  of 
/(a)  will  be  d. 
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Lemma  3.3. 


and 


dim(SW)  <  /( a) 
<  1(a). 


(3.8) 

(3.9) 


Proof 

Fix  a.  To  prove  the  first  part  of  the  lemma  consider  an  arbitrary  7  >  /(a),  and  choose 
rj  >  0  such  that  7  >  /(a)  +  2 rj.  Then,  there  is  £  >  0  and  mo  G  IN  such  that 


iV(n)(a,  e)  <  2n(/(a)+,?) 


for  all  n  >  uiq.  Let  us  define  J(m )  :=  U {kn  :  n  >  m  and  a  —  e  <  .S'j."'1  <  a  +  e}.  Then, 
for  any  m  the  intervals  with  kn  e  J(m)  form  a  cover  of  E^a\  These  intervals  are  of 
length  less  than  2~m.  Moreover,  for  any  m  >  m0  we  have 


| =  Y  N{n\a,e)  ■  2~ni  <  Y  2“n(7_/(a)_??)  <  Y  2”™? 

kn£j(m)  n>m  n>m  n>m 


tends  to  zero  with  m  — >  00.  We  conclude  that  the  7-dimensional  Hausdorff  measure  of 
E [a]  is  zero,  hence,  dimfifb])  <  7.  Letting  7  — >  /(a)  completes  the  first  part. 

Aiming  at  using  /(a)  for  an  estimate  of  Hausdorff  dimensions  consider  an  arbitrary 
7  >  /(a),  and  choose  rj  >  0  such  that  7  >  /(a)  +  2 ?/.  Note  that  A^(r!^(a,e)  <  2n^a?)+r^ 
similar  as  before,  but  this  time  only  for  some  infinitely  many  indices  n  (not  all  large  n) . 
This  is  very  little  information  and  not  sufficient  to  tackle  E^;  even  to  deal  with 
we  need  the  auxiliary  sets  Ki  :=  {t  :  a  —  e  <  Y  <  cl-\-  £  ior  k  =  kn(t)  and  all  n  >  1} 
(this  approach  is  similar  to  [94,  p.  137]).  An  efficient  cover  of  Ki  is  provided  by  the  sets 
I(}"}  where  n  >  l  is  fixed  and  where  k  satisfies  a  —  £  <  s ^  <  a  +  £.  We  find 

Y* '  =  N^n\a,£)  ■  2~ni  <  2 -n<b-Z(a)-r9) 

k :  a— e<s^  <a+e 


where  the  last  inequality  holds  at  least  for  some  infinitely  many  indices  n.  This  allows 
to  conclude  that  the  7-dimensional  Hausdorff  measure  of  Ki  is  zero  and  dirn(Ab)  <  7. 
Since  =  U 1K1  we  conclude  dim(A^al)  <  7  by  countable  continuity  of  the  Hausdorff 
measure  dim(-)  [32,  p.  29].  Letting  7  — ►  /(a)  completes  the  proof.  <0 


3.3  Partition  Function  and  Legendre  Spectrum 

The  second  note  on  the  grain  spectrum  /(a)  concerns  its  interpretation  as  a  Large 
Deviation  Principle  (LDP).  To  this  end  we  consider  N^nl(a,£) /2n  to  be  the  probability 
to  find  (for  a  fixed  realization  of  Y)  a  number  kn  e  {0, . . . ,  2n  —  1}  such  that  Y  G[a  — 
£,  a  +  e].  Typically,  there  will  be  one  “expected”  or  most  frequent  value  of  lim)woos^, 
denoted  o,  and  /(a)  will  reach  its  maximum  1  at  a  =  a.  If  a  differs  from  a,  on  the  other 
hand,  then  [a  —  £,  a  +  e]  will  not  contain  a  for  small  £  and  the  chance  to  observe  coarse 
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exponents  which  lie  in  [a  —  e,a  +  e]  will  decrease  exponentially  fast  with  rate  given 
by  /(a). 

Appealing  to  the  theory  of  LDP-s  we  consider  the  random  variable  An  =  — ns $  hr(2) 
where  K  is  randomly  picked  from  {0, . . . ,  2"  —  1}  with  uniform  distribution  Un  (recall 
that  the  we  study  one  fixed  realization  or  path  of  Y)  and  introduce  its  ‘logarithmic 
moment  generating  function’: 

Definition  3.4.  The  partition  function  of  a  path  of  Y  is  defined  for  all  q  E  JR  as 


r(q)  :=  lirninf  —  —  log2  S(jL\q), 

n— xx)  77, 


where 


2n  — 1 


2"  — 1 


S^n\q)  :=  ^  exp  qns {f]  ln(2))  =  ^  2~nqs^  =  2nJE 


k= 0 


k= 0 


(n) 

2~nqsk 


(3.10) 


(3.11) 


Here,  IEn  stands  for  expectation  with  respect  to  Un.  To  avoid  trivialities  we  set  2  qo°  :=  0 
for  all  q  G  IR,  i.e.,  infinite  s[n')  give  no  contribution  to  the  partition  sum. 

Theorems  on  LDP  such  as  the  one  of  Gartner-Ellis  [27]  apply  then  to  yield  the 
following  result  which  was  established  in  a  slightly  stronger  version  in  [91]: 

Theorem  3.5.  If  the  limit 

r(q)  =  lim  -- log2  S{n\q)  (3.12) 

n— xx)  77, 


exists  and  is  finite  for  all  q  G  IR,  and  if  r(q)  is  a  differentiable  function  of  q, 
double  limit 


f(a)  =  lim  lim  — 

ej.0  n—>oo  n 


log  2N^(a,s) 


then  the 
(3.13) 


exists,  in  particular  f(a )  =  f(a),  and 


/(a)  =  r*(a)  :=  inf  (qa  -  r(q)) 

qeJR. 


(3.14) 


for  all  a. 

For  the  existence  of  the  limit  (3.12)  see  remark  3.11. 

Proof 

The  theorem  of  Gartner-Ellis  [27,  Thm  II]  allows  to  estimate  the  exponential  decay  rate 
of  the  probabilities  Pn[s  f'1  G  E\  of  our  random  variables  s[nl  being  contained  in  a  set  E. 
Recall  that  the  randomness  is  here  in  choosing  the  integer  K  from  {0, . . . ,  2n  —  1}  with 
uniform  distribution  Un.  In  the  light  of  (3.11),  the  assumptions  made  in  the  theorem 
ensure  that  [27,  Thm  II]  is  applicable. 

As  is  typical  for  results  on  large  deviations,  upper  bounds  are  available  for  closed 
sets  E  while  lower  bounds  can  be  obtained  for  open  sets  E.  Recall  that  the  range  of 
allowed  values  for  sf1  in  N(n\a,e)  is  half  open;  thus 

:  —  al  <  <  A^^n)(a,e)  <  ff{k  :  | —  a|  <  e} 
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Applying  the  LDP  bounds  of  [27,  Thm  II]  to  the  above  sets  gives  immediately 

lim  sup  —  log2  (a,  <  sup  r*  (a') 

n—>oo  Tl  \a'—a\<e 

and 

lim  inf  —  log2  (a,  e)  >  sup  r*  ( a ') . 

\a'  —a\<e 

By  continuity  of  r*(a)  these  two  bounds  on  the  right  hand  side  coincide  and  (3.13)  is 
established.  Letting  now  £  — >  0  shows  that  /(a)  =  r*(a).  <0> 

The  assumptions  of  this  theorem  are  too  restrictive  in  many  applications.  Before 
dwelling  more  on  the  relation  between  r  and  /  in  section  4  let  us  note  the  following 
simple  fact.  Its  proof  explains  why  the  Legendre  transform  appears  in  this  context. 

Lemma  3.6.  For  any  a  G  IR 

/(a)  <  r*(a ).  (3.15) 

Proof 

Fix  q  G  IR  and  take  a  with  /(a)  >  —  oo.  Let  7  <  /(a)  and  £  >  0.  Then,  there  are 
arbitrarily  large  n  such  that  N(n\a,  e)  >  2ni .  For  such  n  we  estimate  S^n\q)  by  noting 


2n  — 1 

E(n)  x (n) 

2~nisk  >  V"  2~nqSk 


>  N^(a  £)2~n<'qa+^£')  >  2~n('qa~'y+^£') 


(3.16) 


and  hence  r(q)  <  qa  —  7  +  \q\e.  Letting  £  — >  0  and  7  — »•  /(a),  we  find  r(q)  <  qa  —  /(a). 
Since  this  is  trivial  if  /(a)  =  —  00  we  have  established  that 


r(q)  <  qa  —  f(a )  and  /(a)  <  qa  —  r(q)  for  all  a  and  q  in  IR.  (3-17) 

From  this  it  follows  trivially  that  r(q)  <  f*(q )  and  /(a)  <  r*(a).  <0 

Historical  note  In  the  special  case  when  Y  is  the  distribution  function  At  of  a 
measure  /i,  choosing  the  singularity  exponent  sj”'1  =  ctj^  results  in 

2n  —  l  2"  — 1 

=  E  \M  +  i)2-”)  -  =  27  (f-A”1))"  .  (3.18) 

fc=o  k= 0 

This  is  the  original  form  in  which  r(q)  as  been  introduced  in  multifractal  analysis 
[45,  46,  36,  64],  Note  that  there  is  a  close  connection  to  the  thermo-dynamical  formalism 
[101]. 


3.4  Deterministic  Envelopes 

Often,  we  would  like  to  use  a  simple  analytical  approach  in  order  to  gain  intuition 
on  the  various  spectra  on  a  typical  path  of  Y,  or  at  least  some  estimate  of  them.  To 
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this  end,  we  consider  now  position  on  the  time  axis,  i.e.,  t  or  kn ,  and  the  path  of  Y 
to  be  random  simultaneously  as  we  apply  the  LDP.  More  precisely,  we  consider  the 
exponents  s^\u)  now  as  being  random  variables  over  (Q  x  2n )  where  k  is  picked  with 
uniform  distribution  from  {0, , . . . ,  2”  —  1},  and  independently  of  u. 

Definition  3.7.  The  deterministic  partition  function  of  Y  is 

T(q )  :=  liminf--log2IEn[5'(n)(g)].  (3.19) 

n— xx)  77, 

Remark  3.8.  (Ergodic  Processes)  In  the  definitions  of  r(g)  and  T(q)  we  have 
assumed  that  Y  is  defined  on  a  compact  interval  which  we  took  to  be  [0, 1]  without 
loss  of  generality.  For  processes  defined  on  IR  we  modify  (q)  to 

N2n  —  1 

S{n\q )  :=  lim  —  V  2~nqs 

iV— >00  N  '  ^ 
k= 0 

and  N^n\a,e)  similarly.  For  ergodic  processes  this  becomes  S^n\q)  =  2nWjQ/2~nqsk 
almost  surely.  Thus,  IEn[S'^(g)]  =  S^n\q)  a.s.  and 

T{q)^r{q,u).  (3.20) 

We  refer  to  (7.11)  and  (5.32)  for  an  account  on  the  extent  to  which  marginal  distribu¬ 
tions  may  be  reflected  in  multifractal  spectra  in  general.  4k 

For  processes  on  [0, 1]  we  can  not  expect  to  have  (3.20)  in  all  generality.  Nevertheless, 
(3.20)  holds  in  various  nice  situations  as  we  are  going  to  see,  and  T(q)  does  always  serve 
as  a  deterministic  envelope  of  r(q,u): 

Lemma  3.9.  For  clarity,  we  make  the  randomness  of  r  explicit  by  writing  r(q,u). 
With  probability  one 


r(q,u))  >  T(q)  for  all  q  with  T(q)  <  oo.  (3-21) 

The  inequality  may  be  strict  (see  Example  5.3). 

Proof 

Consider  any  q  with  finite  T(q)  and  let  £  >  0.  Let  n0  be  such  that  IEn[S'^(g)]  < 
2 for  all  n  >  n0'  Then, 


IE 


lim  sup  2nmq)~2£)S{n\q,uj ) 


<  IE  2n{T{q)~2£)S{n\q,co)  <  2~H£  <  00 ■ 

n>no  n>no 


Thus,  almost  surely  lirnsup^^  2n^T^~2e^ S^n\q,  u>)  <  oo,  and  r(q)  >  T(q )  —  2e.  It 
follows  that  this  estimate  holds  with  probability  one  simultaneously  for  all  £  =  1/m 
(' m  G  IN)  and  some  countable,  dense  set  of  q  values  with  T(q)  <  oo.  Since  r(q)  and 
T(q)  are  always  concave  due  to  corollary  4.3  below,  they  are  continuous  on  open  sets 
and  the  claim  follows.  <0> 
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Remark  3.10.  (Importance  of  moments  of  negative  order  q  <  0)  In  ‘tradi¬ 
tional’  statistics  moments  are  usually  considered  to  be  taken  with  respect  to  a  centered 
random  variable,  i.e.,  1E[|X  —  IELA|9].  In  this  setting,  moments  of  negative  order  measure 
only  the  fluctuations  around  the  mean  and  can,  therefore,  often  be  neglected. 

As  we  pointed  out  in  remark  3.8  S^n\q)  can  be  considered  as  approximating  marginal 

( n ) 

moments  of  the  random  variables  2~ns *  ,  at  least  under  some  ergodicity  assumption. 
Depending  on  one’s  choice  these  may  be  moduli  of  increments  as  in  (3.18)  or  moduli  of 
wavelet  coefficients.  Increments  as  well  as  wavelet  coefficients  are  clearly  centered  for 
a  process  with  zero  mean  increments  and  possess  considerable  mass  around  zero,  espe¬ 
cially  for  Gaussian  processes  such  as  ffim;  in  this  case,  negative  order  moments  provide 
indeed  only  little  information.  We  will  comment  in  greater  detail  on  the  multifractal 
scaling  of  fBm  and  the  role  of  negative  q  in  Section  7.3  when  the  necessary  results  are 
available. 

For  processes  with  positive  increments  such  as  cascades,  on  the  other  hand,  negative 
order  moments  become  important  and  relevant,  since  they  capture  the  probability  of 
very  small  increments.  In  other  words,  the  negative  order  moments  are  related  to  the 
time  instances  t  with  high  regularity,  i.e.  the  smooth  parts  in  these  otherwise  ‘spiky’ 
processes. 

A  difficulty  arises  for  cascades  with  fractal  support.  Due  to  boundary  effects  the 
coarse  singularity  exponents  may  become  exceptionally  large  (yet  finite),  causing  the 
partition  function  to  degenerate  for  negative  q  [91].  In  Section  5.6  we  will  show  how 
to  improve  on  the  analysis  using  increment  exponents  afl>  or  Holder  exponents  Iifl> . 
Similar  problems  are  encountered  also  with  wavelet  exponents,  where  a  remedy  has 
been  devised  in  [78]  using  local  maxima  in  wavelet  bands.  4k 

Remark  3.11.  (Quenched  and  annealed  averages) 

A  simple  application  of  Chebichev’s  inequality  shows  that  IE[r(g)]  >  T(q)  which  is 
clearly  not  as  strong  as  lemma  3.9.  However,  IE[r(g)]  is  of  interest  in  itself.  Assuming 
that  the  limit  (3.12)  actually  exists,  Dini’s  theorem  allows  to  exchange  expectation  and 
limit  and  we  may  write 

!E[r(g)]  =  lim  (— l/n)IE[log(S'(n)(g))]. 

n— >  oo 

In  material  science,  this  expression  is  also  known  as  a  quenched  average.  Exchanging 
expectation  and  logarithm  -an  operation  which  in  general  changes  the  object-  we 
obtain  T(q),  also  termed  annealed  average. 

The  free  energy  is  said  to  have  the  self- averaging  property  if  quenched  and  annealed 
averages  are  equal.  Since  r(g)  >  T(g)  almost  surely  the  free  energy  occurring  naturally 
in  the  framework  of  multifractal  analysis  is  self- averaging  if  and  only  if  r(g)  =  T(q) 
almost  surely  (provided  the  limits  exist). 

The  existence  of  the  limit  (3.12)  for  binomial  cascades  has  been  established  in  [18] 
as  well  as  in  [38] .  It  follows  also  from  the  following  simple  observation,  which  promises 
wider  applicability: 

lim  sup  —  —  log2  S^n\q)  <  f*(q ). 

n— xx)  ^ 
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This  fact  can  be  established  similarly  as  lemma  3.6.  Thus,  the  equality  f(a)  =  T*(a ) 
entails  the  existence  of  the  limit  (3.12).  4k 

The  step  from  the  partition  function  r(g)  to  the  deterministic  envelope  T(g)  consists 
in  replacing  averages  of  exponents  within  a  path  by  averages  within  and  across  paths. 
For  /(a)  this  translates  to  replacing  probability  over  one  path  by  the  probability  within 
and  across  paths.  As  we  shall  see  this  means  to  average  A^n)(a,e)  over  all  paths.  For 
simplicity  of  the  argument  fix  n,  a  and  £,  and  let  1  be  the  random  variable  which  is 
1  if  a  —  £  <  (uj)  <  a  +  £  and  0  otherwise,  where  (c o,k)  are  randomly  chosen  from 
0  x  2n.  Obviously,  E[l]  =  P[  1  =  1].  To  compute  this  value  we  make  use  of  Fubini’s 
theorem  and  average  1  first  within  a  fixed  path  -which  yields  N^n\a,  £)  /2n-  and  then 
average  over  all  paths.  Alternatively,  we  may  fix  the  location  k  and  average  1  over  all 
paths  first  -which  yields  Pn[a  —  £  <  (u)  <  a  +  e]-  and  then  average  over  all  k.  In 

summary: 

2n  — 1 

Pnx  2- [a  -.e  <  4n)  <  a  +  e]  =  JEn[Mn\a,£)/2n]  =  2~n  ^  Pn[a  -£<  s <n)  <  a  +  e]. 

k= 0 

(3.22) 

In  analogy  with  (3.5)  we  multiply  this  probability  with  2n  when  defining  the  corre¬ 
sponding  spectrum: 

Definition  3.12.  The  deterministic  grain  spectrum  of  Y  is 

F(a )  :=  lim  lim  sup  —  log2  £)]  (3.23) 

efo  n— »oo  n 

To  have  some  control  over  the  convergence  in  n,  which  is  needed  to  obtain  a  formula 
for  f(a)  valid  for  almost  all  paths  in  Section  8,  we  introduce 

F(a)  lim  lim  inf  —  log2  En[AT^(a,  e)].  (3.24) 

e|0  n— »oo  Tl 

Replacing  Afln^(a,e)  by  (3.22)  in  the  proof  of  theorem  3.5  and  taking  expectations 
in  (3.16)  we  find  properties  analogous  to  the  pathwise  spectra  r  and  /: 

Theorem  3.13.  For  all  a  G  IR 


F(a)  <  T*(a).  (3.25) 

Furthermore,  ifT(q)  admits  a  finite  limit  as  n  — >  oo  for  all  q  G  IR  similar  to  (3.12), 
and  is  concave  and  differentiable  as  a  function  of  q,  then  F(a)  admits  a  limit  as  n  oo 
analogous  to  (3.13),  in  particular 

F(a)  =  F(a)  =  T*(a).  (3.26) 

We  give  such  an  example  in  Example  5.3. 

It  follows  from  lemma  3.9  that  with  probability  one  T*(a,uj)  <  T*(a)  for  all  a. 
Similarly,  the  deterministic  grain  spectrum  F(a)  is  an  upper  bound  to  its  pathwise 
defined  random  counterpart  /(a, uj),  however,  only  pointwise.  On  the  other  hand,  we 
have  here  almost  sure  equality  under  certain  conditions. 
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Theorem  3.14.  Fix  some  number  a  G  IR.  Then,  almost  surely 

f(a,  u>)  <  F(a).  (3.27) 

Ifforalln  the  random  variables  (k  =  0, . . .  ,2n  —  1)  arei.i.d.,  andifF(a)  =  F(a)  > 
0,  then  almost  surely 

/(a,  u)  =  /(a,  u)  =  F(a).  (3.28) 

Compare  the  regularity  condition  F(a)  =  F(a)  (see  (3.26))  to  the  more  restrictive 
requirement  that  F(a)  assumes  a  limit  similar  to  (3.13). 

Remark  3.15.  (Independent  increments)  It  is  easy  to  extend  this  result  and 
allow  s ^  to  depend  on  some  of  its  nearest  neighbors,  say  on  sfl>  for  \l  —  k\  <  m0  for 
some  constant  mo-  Thus,  if  Y  has  independent  increments,  (3.28)  applies  not  only  to 
the  increment  exponents  but  also  to  the  Holder  exponents  /q”'1  as  well  as  to  the 
wavelet  exponents  wjp  for  compactly  supported  wavelets.  4k 

Proof 

The  inequality  (3.27)  follows  as  in  lemma  3.9,  using  the  estimate 

IEnlimsup2“n(F(a)+2£)iV(n)(a,e)  <  ^  2-n(F(a)+2e)iV(n)(a,  e)  <  2“n£. 

n— >oo  .  . 

n>m0  n>m  o 

Since  the  grain  spectra  are  not  necessarily  continuous,  the  inequality  cannot  be  estab¬ 
lished  for  all  a  simultaneously,  but  only  for  a  countable  set  of  a-values. 

Assume  now  that  the  .s[n'1  [k  —  0, . . . ,  2n  —  1)  are  independent  and  identically  dis¬ 
tributed.  To  show  equality  in  (3.27),  let  us  note  first  that  N^n\a,£)  is  a  Bernoulli 
variable: 

Pn  [N™(a,  e)  =  j]  =  (  T.  )  pj(  1  -  pn)2^  (3.29) 

where  pn  =  Pn[a  -  £  <  s <  a  +  e]  =  2_nlE^[A^(n)(a,  e)]. 

The  property  (3.26)  says  that  the  Pn[a  —  £  <  s^';>  <  a  +  e]  are  ‘close  to’  converging  as 
(n  — *  oo).  More  precisely,  (3.26)  guarantees  that  for  any  p  >  0  we  find  £o(p)  such  that 
lirn  supn  and  liminfn  of  these  quantities  do  not  differ  by  more  than  p  for  all  £  <  £0. 
Thus,  for  any  such  £  and  any  rj  >  0  there  is  711(77,  s,  p)  such  that  for  all  n  >  ri\ 

F(a)  +  p  +  p  >  -  log2(2 npn)  >  F(a)  -77  -  p.  (3.30) 

n 

Let  now  p  >  0  and  77  >  0  be  such  that  1  >  F(a)  +  77  >  F(a)  —  77  >  0.  Using  (3.30) 
it  follows  easily  that  Pq  [iV(”)(a,  e)  =  j]  / Pq  [AT^(a,£)  =  j  —  l]  >1,  i.e.,  (3.29)  grows 
monotonously  as  a  function  of  j,  for  j  <  2nF(ah 

Now,  let  0  <  7  <  F(a)  —  77  —  p  and  choose  l  such  that  l  —  1  <  2n7  <  l.  Then, 
exploiting  the  monotony  of  (3.29),  one  finds 

Pn  [N^n\a,£)  <  2n7]  <  l  •  (  2"  )  pn\l  -  pn)2n~l .  (3.31) 
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Using  a  standard  estimate  on  binomial  coefficients  based  on  Stirling’s  formula  and 
observing  that  F(a)  —  rj  —  p  >  0a  tedious  but  straightforward  calculation  allows  to 
conclude  that  the  RHS  of  (3.31)  decays  with  hyper-exponential  speed  in  n,  mainly 
because  of  the  last  term.  In  summary, 

OO 

^Pn  [N^{a,e)<  2n7]  <  oo. 

n=  1 

By  the  Borel-Cantelli  lemma 

Pci  [N^n\a,  e)  <  2n7  for  infinitely  many  n]  =  0.  (3.32) 

In  other  words,  almost  surely  N^n\a,e)  >  2n7  for  all  large  n,  or  again  in  other  words, 
almost  surely  liminfn_>00(l/n)  log2  N^n\a,  e)  >  7.  Letting  £  — >  0  gives  f(a,u>)  >  7  a.s. 
Letting  then  7  — >  F(a)  —  rj  —  p,  rj  —>■  0  and  p  — >  0  (all  in  discrete  sequences,  of  course) 
yields  that  almost  surely 

F(a)  <  f(a,u). 

But  /(a,  a;)  <  f(a,u)  by  definition  and  f(a,u)  <  F(a)  a.s.  by  (3.27).  So,  (3.28)  follows. 

❖ 

Here  is  a  first  result  that  allows  to  compute  almost  sure  pathwise  spectra  from 
knowing  only  T  under  various  regularity  assumptions. 

Corollary  3.16.  Assume  that  T(q)  admits  a  finite  limit  as  n  — >  00  for  all  q  G  1R,  and 
is  concave  and  differentiable  as  a  function  of  q.  Assume  furthermore  that  for  large  n 
the  singularity  exponents  s^l>  (k  —  0, .  2n  —  1)  used  in  T  are  i.i.d.  random  variables. 
Pick  a  such  that  T*(a)  >  0.  Then  almost  surely 

f(a,u)  =  f(a,u)  =  T*(a,cu)  =  F(a )  =  T*(a).  (3.33) 

Remark  3.17.  (Negative  Dimensions) 

Note  that  T*  and  F  may  assume  negative  values,  which  is  not  possible  for  /.  Con¬ 
sequently,  T*  and  F  may  be  expected  to  be  a  good  estimator  of  /  only  where  they  are 
positive. 

Negative  F(a)  and  T*(a )  have  been  termed  negative  dimensions  [67].  They  corre¬ 
spond  to  probabilities  of  observing  a  coarse  Holder  exponent  a  which  decay  faster  than 
the  2"  ‘samples’  of  available  in  one  realization.  Oversampling  the  process,  i.e.,  an¬ 
alyzing  several  independent  realizations  will  increase  the  number  of  samples  and  more 
‘rare’  s'f  '>  may  be  observed.  In  loose  terms,  in  exp(— nln(2)F(a))  independent  traces 
one  has  a  fair  chance  to  see  at  least  one  sfi'1  of  size  ~  a. 

Thereby,  it  is  essential  not  to  average  the  spectra  /(a)  of  the  various  realizations  but 
the  numbers  N^n\a,e).  This  way,  negative  ‘dimensions’  /(a)  become  visible.  4k 


4  The  Multifractal  Formalism 

Various  multifractal  spectra  have  been  introduced  in  the  previous  section,  along  with 
some  simple  relations  between  them  which  we  may  summarize  as: 
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Corollary  4.1  (Multifractal  formalism). 

Provided  all  spectra  are  in  terms  of  the  same  singularity  exponent,  we  have 

clim(A^)  <  f(a)  <  /(a)  <  r*(a)  <'  T*(a).  (4.1) 

The  first  relations  hold  pathwise,  and  the  last  one  with  probability  one  for  every  a  G  IR. 
Similarly, 

clim(£[o])  <  /(a)  <’  F(a)  <  T*(a). 

The  spectra  on  the  left  end  have  stronger  implications  on  the  local  scaling  structure 
while  the  ones  on  the  right  end  are  more  easy  to  estimate  or  calculate. 

This  set  of  inequalities  could  fairly  be  called  the  ‘multifractal  formalism’.  However, 
in  the  mathematical  community  a  slightly  different  terminology  is  already  established 
which  goes  as  ‘the  multifractal  formalism  holds’  and  means  that  for  almost  all  paths 
of  a  particular  process  the  spectrum  dim(A^al)  coincides  with  the  Legendre  transform 
of  some  adequate  partition  function  (such  as  r(q)).  It  appears  that  this  property  holds 
indeed  in  a  generic  sense,  meaning  that  in  the  proper  context  (see  [52])  the  ‘multifractal 
formalism’  is  valid  quasi-surely  in  Baire’s  sense.  In  view  of  (4.1)  this  implies  that  equality 
holds  then  between  all  introduced  spectra. 

Though  we  pointed  out  some  conditions  for  equality  between  /,  t*  and  T*  we  must 
note  that  in  general  we  may  have  strict  inequality  in  some  or  all  parts  of  (4.1).  Such 
cases  have  been  presented  in  [91]  and  [94],  There  is,  however,  one  equality  which  holds 
always  and  connects  the  two  spectra  in  the  center  of  (4.1). 

Theorem  4.2.  Consider  a  realization  or  path  ofY .  Recall  that  infinite  s'^}  don’t  con¬ 
tribute  to  r(q)  nor  to  f(a). 

a)  Both-sided  multifractal:  If  the  finite  sfi'1  are  bounded,  then 

r(q)  —  /*(<?)  for  Q  £  IT  (4.2) 


b) 


Left-sided  multifractal: 

from  below,  then 


If  the  finite  s 1"'> 


r(g)  = 


/*(?) 

— oo 


are  unbounded  from  above  but  bounded 

for  all  q  >  0 
for  all  q  <  0. 


c)  Right-sided  multifractal:  If  the  finite  sfi'1  are  bounded  from  above  but  unbounded 
from  below,  then 

—oo  for  all  q  >  0 
f*(q)  for  all  q  <  0. 


r(g)  = 


d)  If  the  finite  sj."'1  are  unbounded  from  above  and  from  below,  then 


r(q)  =  —oo 


for  all  q  0. 
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Proof 

The  following  notation  will  be  useful:  Kn(a,e )  :=  {k  :  a  —  £  <  ,s[n')  <  a  +  e}.  Recall 
that  r(q)  <  f*(q)  from  lemma  3.6. 

Now,  to  estimate  r(q)  from  below,  we  will  group  the  terms  in  S^n\q)  conveniently, 
i.e., 

\a/e\ 

E  £  +  £ 

i=-[a/e\  Kn(ie,e)  |s<>)|>s 

where  we  keep  the  choice  of  a  open  for  the  moment.  Since  we  need  uniform  estimates 
on  j\dn)  (a,  e)  for  various  a,  some  preparation  is  needed. 

Fix  q  G  IR  and  let  r/  >  0.  Then,  for  every  a  G  [—a,  a]  there  is  £o(a)  and  no  (a)  such  that 
N^(a,s)  <  2"CO)+»?)  for  all  £  <  £o(a)  and  all  n  >  n0(a).  We  would  like  to  have  e0  and 
no  independent  from  a  for  our  uniform  estimate.  To  this  end  note  that  N^n\a',£')  < 
N^n\a,£)  for  all  a'  G  [a  —  £/ 2,  o  +  e/2]  and  all  £r  <  e/2.  By  compactness  we  may  choose 
a  finite  set  of  a,j  ( j  =  1, . . . ,  m)  such  that  the  collection  [aj  —  £o(a.j)/2,  %  +  £o(«j)/2] 
covers  [—a,  a].  Set  £\  =  (1/2)  minj=ir..jm  £oiflj)  and  ri\  =  maxj=ir..)rn  no(a,j).  Then,  for 
all  £  <  £i  and  n  >  ri\  we  have 


(n) 

-nqs\ 


(4.3) 


S{n\q)  <  ( 


N(n\a,£)  <  2n(/(a)+,?) 


for  all  a  G  [—a,  a].  For  £  <  £i  and  n  >  n±  we  estimate  the  first  term  in  (4.3)  by 


La/eJ  \a/e\ 

2~nqs<k)  <  N^n\i£,£) 2-n(9*£”bk)  < 

i=—  \a/s\  Kn (ie,e)  i=—\a/s J 

\a/e\ 

2_n(9*£-/(*£)-J?-kle)  <  (2|_a/£j  +  1)  •  2~n^*^~v~^£\  (4.4) 

i=—  [a/e  J 

Case  a):  For  bounded  s ^  this  is  all  we  need.  Indeed,  choosing  a  larger  than  s[n')  for 
all  n  and  k,  the  second  term  in  (4.3)  vanishes  and  (4.4)  estimates  S(n\q)  itself.  Letting 
n  — >  oo  we  find  r(q)  >  f*(q )  —  q  —  \q\£  for  all  £  <  £i.  Now  we  let  £  — >  0  and  finally 
r)  — >  0  to  find  the  desired  inequality  r(q)  >  f*(q ). 

Case  b):  If  the  are  bounded  only  from  below  we  proceed  differently  for  positive 
and  negative  q.  For  q  >  0  we  choose  a  large  enough  to  ensure  qa  >  f*(q)  +  1  as  well  as 
,s[n)  >  —a  for  all  k  and  n.  The  second  term  in  (4.3)  is  then  bounded  by 

2~nqs(£)  <  2n2~nqa  <  2~nk*^ 


using  that  the  sum  has  at  the  most  2n  terms.  This  expression  is  certainly  smaller  than 
the  right  hand  side  of  (4.4),  whence  S^n\q)  is  bounded  by  twice  (4.4)  and  the  result 
follows  as  before  in  case  a). 


4  The  Multifractal  Formalism  33 


For  q  <  0  we  simply  note  that  for  any  number  x  we  can  find  arbitrarily  large  n  such 
that  Sfc  >  x  for  some  k.  This  implies  that  S('n\q)  >  2~nqx  and  r(q)  <  qx.  Letting 
x  — >  oo  proves  the  claim  r(q)  =  — oo  since  q  is  negative. 

Case  c):  If  the  s/l>  are  bounded  only  from  above  we  argue  similarly  as  in  case  b). 
Finally,  in  case  d)  we  find  r(q)  =  — oo  for  all  q  fz  0  as  in  case  b)  and  c).  <0> 


Note  that  r(q)  may  be  discontinuous  at  0  in  the  one-sided  cases,  but  as  a  Legendre 
transform^  it  is  always  concave.  Moreover: 

Corollary  4.3  (Properties  of  the  partition  function).  The  partition  function  r(q) 
is  always  concave.  Depending  on  which  case  of  theorem  f.2  applies,  r(q)  is  moreover 
continuous  on  1R,  {q  >  0},  or  {q  <  0};  furthermore,  it  is  differentiable  in  those  sets 
with  at  the  most  countable  many  exceptions. 

If  for  large  enough  n  all  s are  positive  for  a  given  path  of  the  process  Y,  then  r(q) 
is  non- decreasing.  Such  is  the  case  for  s^  =  a ^  and  s//'>  =  h provided  the  path  is 
increasing,  or  at  least  of  bounded  variation,  and  also  for  =  w ^  provided  that  the 
path  lies  in  an  appropriate  regularity  space  such  as  L2,  the  space  of  square- integrable 
functions.  However,  r(q)  may  have  decreasing  parts  if  the  analyzed  process  exists  only 
in  the  distributional  sense,  such  as  binomial  measures  (see  section  5  and  Figure  7). 

From  a  numerical  point  of  view  r(q)  is  more  robust  than  /  since  r(q)  involves  averages 
and  is  not  in  terms  of  a  double  limit.  Thus,  one  would  like  to  invert  theorem  4.2.  In 
order  to  efficiently  do  so  we  need: 

Lemma  4.4  (Lower  semi-continuity  of  f  and  F ).  Let  am  converge  to  a *.  Then 


/(a*)  >  limsup  f(am) 

m— »oo 


(4.5) 


and  verbatim  for  F. 


Proof 

For  all  £  >  0  one  can  find  m0  such  that  a*  —  £  <  am  —  ej 2  <  am  +  e/2  <  a*  +  £  for 
all  m  >  m0.  Then,  N(n\a*,e)  >  N{-n\am,£/2 )  and  !E[lV(n)(a*, £)]  >  IE [lV^(am, e/2)\. 
Following  now  only  the  case  for  /,  we  write 


lim  sup  —  log2  N('nl  (a*,  e) 

n—>  oo  ^ 


>  lim  sup  -  log2  N(n)(am,  e/2 )  >  f(am) 


n 


for  any  m  >  m0(e).  Letting  first  m  — >  oo  and  then  £  — >  0  proves  the  claim.  <) 


Corollary  4.5  (Central  multifractal  formalism).  We  always  have 

/(a)  <  f**(a)  =  r*(a).  (4.6) 

For  some  a  we  have  even  equality.  More  precisely,  for  all  q  G  1R  for  which  r'(g±)  is 
meaningful 

/(a)  =  r*(a)  =  qr\q±)  -  r(q±)  ata  =  r'(q±).  (4.7) 


§For  a  tutorial  on  the  Legendre  transform  see  [89,  App.  A]. 


34 


R.  H.  Riedi,  Multifractal  Processes 


Note  that  (4.7)  uses  only  local  properties  of  r(q)  which  should  be  compared  to 
theorem  3.5. 

For  a  concave,  real- valued  function  g(q)  the  derivative  g'  exists  in  all  but  countable 
many  exceptional  points  and  it  is  a  monotone  function.  Thus,  one  can  define  mean¬ 
ingfully  g'(q+)  for  all  q  as  the  limit  of  g'(qn)  where  qn  approaches  q  from  the  right 
and  runs  through  any  set  of  points  where  the  derivative  exists.  The  only  difficulty  with 
r(q)  lies  then  in  its  possibly  infinite  values.  According  to  theorem  4.2  and  corollary  4.3 
r'(g±)  is  meaningful  either  for  all  q  G  IR,  all  q  >  0,  all  q  <  0,  or  none. 

Remark  4.6.  (Special  values  of  r  and  the  shape  of  /)  By  convention,  S'^(O) 
counts  how  many  of  the  sf1'1  are  finite.  So,  S'(n')(0)  >  0  and  we  have  always  r( 0)  <  0. 
Moreover,  S'^(O)  >  N^n\a,e)  for  all  a,  which  implies 

/(a)  <  — t(0)  (4.8) 

with  equality  at  a  =  r'(0±),  lest  the  finite  sf'}  are  unbounded  from  below  and  above. 
If  all  4n)  are  finite,  then  r(0)  =  —1.  In  general,  —  r(0)  gives  the  (box)  dimension  of  the 
support  of  the  analyzed  path  of  the  process  Y  (compare  Example  5.4  and  [91]). 

In  some  cases,  such  as  sf'1  =  tiff'1 ,  the  singularity  exponents  are  monotonous  in  the 
sense  that  2~nsk }  >  max(2_ns2^\  2_ns2fc+i).  If  so,  S^n\  1)  >  2S'(^n+1^(l)  and  r(l)  >  —1 
( iff  '1  are  always  bounded  from  below).  In  the  ‘classical  case’  of  sf1^  =  otff*  for  an 
increasing  path  of  Y  we  have  obviously  S'fifi(l)  =  S'(n+1)(l)  and  r(l)  =  0.  Consequently, 

Increasing  path:  /«(o)  <  a  (4.9) 

with  equality  at  both,  a  =  r'(l+)  and  a  =  r'(l— )  (see  Figure  6).  It  cannot  be  stressed 
enough  that  (4.9)  may  fail  for  processes  which  are  not  monotonous  (compare  Sec¬ 
tion  7.3).  For  certain  increasing  processes,  r'(  1)  has  been  identified  as  the  dimension  of 
the  carrier  and  received  the  name  ‘information  dimension’  [44,  42,  43,  83,  32,  91,  6,  37] 
(compare  (5.27)).  One  can  think  of  the  support  as  a  function’s  ‘realm’,  i.e  the  closure  of 
all  points  where  the  function  is  non-constant,  and  of  the  carrier  as  the  function’s  ‘core’, 
i.e.,  the  smallest  set  of  points  one  can  choose  over  which  the  function  still  exhibits 
its  full  variation.  For  a  cumulative  distribution  function  of  a  probability  measure,  the 
support  collects  all  but  the  open  intervals  with  zero  probability,  while  the  core  marks 
the  smallest  set  with  probability  one. 

* 

Proof 

The  graph  of  /**  is  the  concave  hull  of  the  graph  of  /  which  implies  (4.6).  It  is  an 
easy  task  to  derive  (4.7)  under  assumptions  suitable  to  make  the  tools  of  calculus 
available  such  as  continuous  second  derivatives.  The  reader  not  familiar  with  Legendre 
transforms  should  take  the  opportunity  and  perform  a  simple  minimization  of  qa  —  r(q) 
as  a  function  of  q  to  convince  himself  of  (4.7). 

To  prove  (4.7)  in  general  let  us  first  assume  that  r  is  differentiable  at  a  fixed  g*.  Since 
r(g*)  —  f*(q*)  we  can  find  a  sequence  am  such  that  r(g*)  =  linim^oo  g*am  —  f(am). 
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Now  consider  the  functions  q  i— >•  qam  —  f(am).  Due  to  (3.17)  they  must  all  lie  above 
t(-).  But  evaluated  at  g*  these  functions  approach  r(g*).  Being  linear,  their  slopes  must 
converge  to  the  slope  of  the  tangent  at  r(q )  in  g*,  in  other  words,  am  must  converge  to 
a*  :=  r'(g*). 

From  this,  f(am )  converges  to  g*a*  —  r(g*).  Applying  lemma  4.4  we  find  that  /(a*)  > 
limm/(am)  =  g*a*  —  r(g*).  Recalling  (3.17)  again  implies  the  desired  equality  /(a*)  = 
g*a*  -  r(g*). 

Let  us  now  consider  any  g  G  IR  for  which  r'(q+)  is  meaningful.  Then,  there  is  a 
sequence  of  numbers  qm  larger  than  q  in  which  r  is  differentiable  and  which  converges 
down  to  q.  Consequently,  r'(q+)  =  lim mr'(gm).  The  formula  (4.7)  being  established 
at  all  qrn  lemma  4.4  applies  now  with  am  =  r'(gm)  and  with  a*  =  r'(g+)  to  yield 
/(r'(g+))  =  /(a*)  >  lirnsupm/(am)  =  limm  gmr'(gm+)  -  r(gm+)  =  gr'(g+)  -  r(g+). 
Again,  (3.17)  furnishes  the  opposite  inequality.  A  similar  argument  applies  to  r'(g— ). 

0 

Corollary  4.7.  If  T(q)  is  finite  for  some  positive  q  then  the  finite  s ^  are  bounded 
from  below  for  almost  all  paths;  moreover,  for  all  q  >  0  with  finite  T(g) 

T(g)  =  F*(q).  (4.10) 

If  T(g)  is  finite  for  some  negative  q  then  the  finite  are  bounded  from  above  for 
almost  all  paths,  and  (4-10)  holds  for  all  q  <  0  with  finite  T(g).  Moreover,  for  all  q  for 
which  T'(g)  exists 


F(a)  =  T*(a)  =  qT'(q)  —  T(g)  ata  =  T\q).  (4.11) 

See  Example  5.3  for  a  case  with  F(a)  <  T*(a). 

Proof 

If  T(g)  is  finite  then  r(g)  >  —  oo  by  lemma  3.9.  So,  the  first  claims  follow  from  theo¬ 
rem  4.2  and  by  taking  expectation  in  (4.3)  and  (4.4).  Formula  (4.11)  could  actually  be 
strengthened  to  read  similar  to  (4.7),  the  proof  following  along  the  lines  of  corollary  4.5. 

❖ 


5  Binomial  Multifractals 

The  binomial  cascade  has  a  long  standing  tradition  in  serving  as  the  paradigm  of 
multifractal  scaling  [64,  56,  66,  15,  49,  10,  91,  95].  We  present  it  here  with  an  eye  on 
possible  generalizations  of  use  in  modelling  increasing  random  processes. 

5.1  Construction  of  conservative  cascades 

The  binomial  cascades  form  a  class  of  increasing  processes.  The  construction  of  a 
binomial  cascade  relies  on  two  main  ingredients:  the  first  one  provides  infinite  geometric 
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detail  through  iterative  sub-division  while  the  second  one  provides  randomness.  To 
clearly  distinguish  these  two  parts  we  construct  a  binomial  cascade  pathwise,  i.e.,  we 
first  describe  the  geometrical  detail  inherent  to  any  path  (or  realization)  of  the  cascade 
and  introduce  randomness  only  afterwards. 

Each  path  of  a  binomial  cascade  is  a  distribution  function  A4b(t),  in  particular  each 
path  is  right-continuous  and  increasing.  As  we  will  see,  almost  all  paths  possess  no 
derivative  (or  density  function),  meaning  that  M.b{t)  f_  Ad(,(s)ds;  unless  for  a 
handful  of  trivial  cascades.  Being  distribution  functions,  though,  each  path  TVT is 
uniquely  related  to  a  measure  or  (probability)  distribution  fib  through  M.b{t)  =  R&(]  — 
oo,f[).  Actually,  the  easiest  way  to  construct  M.b  is  by  defining  pb  first.  Except  for 
trivial  cases,  pb  is  a  true  distribution  and  not  a  function  in  the  usual  sense.  We  will 
call  /if,  a  binomial  measure  and  M.b{t)  =  pb(]  —  oo,£[)  a  binomial  cascade. 

In  order  to  define  pb  we  make  the  notation  (2.4)  more  precise  in  the  following. 
For  any  fixed  t  there  is  a  unique  sequence  ki,k2,...  such  that  the  dyadic  intervals 
Ikn  =  [kn2~n,  ( kn  +  l)2~n[  contain  t  for  all  integer  n,  indeed,  kn  =  kn(t )  =  |_£2nJ.  It  is 
clear  that  any  t  defines  in  this  way  a  unique  sequence  kn  (n  G  IN)  with  the  following 
property: 

Definition  5.1.  We  call  a  sequence  (£q,  &2,  •  •  •)  such  that  is  a  subinterval  of  1^ 

a  nested  sequence. 

Vice  versa,  for  an  infinite  nested  sequence,  the  /["'*  form  a  decreasing  sequence  of 
half  open  intervals  which  shrink  down  to  a  unique  singleton  {£}. 

The  idea  behind  a  binomial  measure  is  to  redistribute  the  mass  lying  in  ij^  among 
its  two  dyadic  subintervals  l!f^ 1  and  I^++  i  in  the  proportions  of  certain  given  numbers 
and  M^^\.  For  consistency  we  require  M^1'*  +  =  1.  To  make  this 

procedure  meaningful  and  amendable  to  analysis,  further  restrictions  will  be  imposed 
on  the  multipliers  M^l>  in  an  instant  (see  (i)-(iii)  below). 

Definition  5.2.  For  given  multipliers  M^"'1  the  binomial  measure  is  defined  by  setting 


w(4:>)  =  ■  M/y1’  ■  ■  ■  m<;>  ■  m‘0),  cs.i) 

and  the  binomial  cascade  is  given  by  setting  A4&(0)  =  0  and 

Mb((kn  +  1)2“")  -  Mb{kn2~n)  =  pb(I^).  (5.2) 


It  is  enough  to  define  the  mass  (or  probability)  of  dyadic  intervals  since  any  interval 
]  —  oo,£[  can  be  written  as  a  disjoint  union  of  dyadic  intervals  J ^  and  A4b(t)  = 
pb(]  —oo,t[=  Y2nfJ,b(J^)-  particular,  integrals  (expectations)  with  respect  to  fib  can 
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FIGURE  5.  Iterative  construction  of  the  binomial  cascade. 


To  say  the  same  in  other  words,  (5.2)  defines  Mb  in  all  dyadic  points.  As  a  dis¬ 
tribution  function,  Mb  is  continuous  from  the  right  side.  Thus,  knowing  the  function 
Mb  in  dyadic  points  is  sufficient.  Note  that  Mb  is  continuous  at  a  given  point  t  unless 

=  1  f°r  all  n  large. 

To  generate  randomness  in  Mb,  we  choose  the  various  to  be  random  variables. 
The  above  properties  hold  then  pathwise.  We  will  make  the  following  assumptions  on 
the  distributions  of  the  multipliers  M 

(o):  Total  mass.  Any  suitable,  positive  random  variable  may  be  chosen  as  Mq°\ 

(i) :  Conservation  of  mass.  Almost  surely  M^l>  is  positive  and 

M{£n  +  M'tll  =  1  (5.5) 

for  all  n  and  k.  As  we  have  seen,  this  guarantees  that  Mb  is  well  defined. 

(ii) :  Nested  independence.  All  multipliers  of  a  nested  sequence  are  mutually  inde¬ 

pendent.  As  a  consequence,  for  any  nested  sequence  as  in  (5.2)  we  have 

■  ■  ■  M<°>]  =  En[M<;>] . . .  En[M<0)]  (5.6) 

and  similar  for  other  moments.  This  will  allow  for  simple  calculations  in  the 
sequel. 
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(iii):  Identical  distributions  For  all  n  and  k 

A  more  general  version  of  (iii)  was  given  in  [89]  to  allow  for  more  flexibility  in  model 
fitting  (see  Example  7.2). 

5.2  Relaxed  conservation 

The  theory  of  cascades  or,  more  properly,  T-martingales  [56,  10,  49,  8],  provides  a 
wealth  of  possible  generalizations.  Most  importantly,  it  allows  to  remove  the  depen¬ 
dence  between  sibling  multipliers  enforced  by  the  almost  sure  conservation  condition 
(5.5)  and  to  require  only  that 

(i’):  Conservation  in  the  mean 


EA  [Mg  +  Mi]  —  1. 


(5.8) 


Relaxing  (i)  into  (i’)  brings  about  various  peculiarities.  Some  of  these,  such  as  virtual 
scaling  exponents  require  an  appropriate  context  to  be  fully  appreciated  and  we  defer 
their  discussion  to  Example  5.3.  Here,  we  relate  merely  the  changes  in  the  construction 
that  are  required.  Indeed,  since  the  total  increment  is  not  preserved  under  iteration 
the  definition  of  the  cascade  (5.1)  needs  to  be  adapted  to  read  as: 


milP)  =  K’ ■  Mi 


(n) 

kn 


jr(n- 1) 
kn  —  l 


m<;'  .  Mf 


lim  V"  M, 

m. — ky)  •  ^ 


(n+l) 

2kn-\-ii 


M, 


(n+2) 


4fcn+»2 


M. 


(n+m) 


2rnkn+i.n 


(5.9) 


Here,  the  sum  runs  over  all  nested  sequences  (H, . . . ,  im)  of  length  m.  This  sum  reduces 
to  1  in  the  case  of  strict  conservation  (5.5),  which  follows  by  induction  in  m.  Similarly, 
the  expectation  of  the  inner  sum  equals  1  due  to  (i’).  Moreover,  the  sum  forms  a  positive 
martingale  and  as  such  converges.  This  makes  (5.9)  meaningful. 

The  limit,  however,  may  degenerate  and  vanish  almost  surely.  This  may  happen 
despite  the  conservation  of  mass  in  the  mean,  because  products  of  multipliers  might  be 
very  small  with  quite  large  probability  and  the  binomial  cascades  ‘dies  out’.  This  case 
is  equivalent  with  IE[.M&(1)]  =  0,  since  A4b(l)  is  a  positive  random  variable,  as  well  as 
equivalent  with  T'(l)  <  0  (see  [56]).  For  an  intuitive  reasoning  we  mention  that  T'{  1) 
is  the  almost  sure  dimension  of  the  carrier  of  the  cascade  whenever  it  is  positive  (see 
[56],  remark  5.8  and  Example  5.3).  Let  us  assume  for  the  remainder  that  T'{  1)  >  0; 
then  IE[A4b(l)]  =  1E[Mq0)]. 

We  noted  that  Aib(  1)  is  in  general  not  equal  to  Mq°\  Similarly,  /^(J^)  is  not  equal 
to  -M(S0)  in  general.  The  ‘identically  distributed  multipliers’  (ii),  however, 

imply  that  Hb(lj^)  is  distributed  as 


(5.10) 
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Mandelbrot  [66]  calls  the  low  frequency  part  of  pb(ljf^),  and  the  subtree 

of  multipliers  ‘hanging  off’  below  the  node  kn  -which  essentially  constitute  the  limiting 
part  distributed  as  M.b(l)~  the  high  frequency  part.  Since  !E[A/f6(l)]  =  E[Mo°^  we  get 

e[»(4:')  i  «£>, , , , ,  m<;>]  =  m<”>  ■  ■  ■  m<;>  .  e[m<°>].  (s.hj 

and  Elftf/';1)]  =  1/2"E[M<01]. 

For  numerical  simulation  (5.11)  reminds  us  that  an  iterative  computation  of  the  prod¬ 
ucts  does  not  provide  the  value  of  the  actual  process  increment  pb(l[fL)) 

-as  is  the  case  with  strict  conservation  (5.5)-  but  rather  what  its  value  could  be  expected 
to  be  if  the  construction  continued. 

The  main  advantage  of  relaxing  (5.5)  is  that  we  can  use  unbounded  multipliers  Mo 
and  Mi  such  as  log-normal  random  variables.  In  this  particular  case,  the  marginals  of 
the  increment  process,  i.e.,  pb(ljfi^)  are  exactly  log-normal  on  all  scales.  For  general 
binomials  with  (i’),  (ii)  and  (iii)  it  can  be  argued  that  the  marginals  pb(ljfi^)  are  at 
least  asymptotically  log-normal  by  applying  the  Central  Limit  Theorem  provided  the 
logarithm  of  (5.2)  is  of  finite  variance. 


5.3  Wavelet  Decomposition 


The  scaling  coefficients  of  pb  using  the  Haar  wavelet  are  simply 


Haar:  Dn,k(p b)  =  /  <f*)k[t)  pb(dt)  =  2n/2 


»(fc+i)2- 


'  k2~ 


W,(dt)  =  2”/V(4”1)  (5.12) 


from  (2.8)  and  (5.3).  With  (2.9)  and  (5.2)  we  get  the  explicit  expression  for  the  Haar 
wavelet  coefficients: 


Haar: 


2-n/2Cn,„M  =  iML1’)  -  wUSi)  = 


id) ( i\/r(n+^)  _  /i/fC+b  a 
2 kn  lvl2kn+l) ' 


i= 0 


(5.13) 

To  obtain  formulas  for  general  wavelets  let  us  start  with  the  observation  that  for  any 
function  if  supported  on  [0, 1] 

2 -n/2C„, *„(«.)=  /,  y  (2<">t  -  W»(d()  =  M<”> . . .  M<‘> .  f  (5.14) 


/■(«) 

lkn 


This  comes  about  from  (5.2)  and  (5.3).  Here  p\fl'kn^  is  the  binomial  measure  constructed 
from  the  subtree  which  has  as  its  root  at  the  node  kn  of  level  n  of  the  original  tree: 


Definition  5.3.  A  binomial  sub-cascade  is  defined  for  any  given  nested  sequence 
ki, ...  ,kn  by  setting 


kb 


(4m))  =  M0(0)M 


(n+1) 
2&n_Ml 


j\An+ 2)  i\An+m) 

'  1V1Akn+l2  '  '  '  1V12 mkn+in 


(5.15) 


for  any  nested  sequence  (H,  •  •  •  Am),  where  Mq0'1  is  an  independent  copy  of  Mq0) .  Also, 
we  set  M.bl'kn\t)  :=  p^’kri\[0,t]).  The  generalization  to  cascades  with  conservation  in 
the  mean  is  obvious,  but  tedious. 
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With  little  more  effort  we  compute  the  analogue  to  (5.14)  for  Mb  itself.  In  preparation 
let  us  note  that  for  kn  =  kn(t)  =  [t2n\ 

Mb{t )  -  Mb{kn2~n)  =  iib([kn2~n,t})  =  M%>  ■  ■■M^Mtl’kn)(2nt  -  kn).  (5.16) 

Let  us  assume  now  that  f  ip  =  0  and  that  the  function  ip  is  supported  on  [0, 1],  whence 
ip2n  •  —k)  is  supported  on  I^l> .  Substituting  t'  =  2nt  —  kn  yields 


J  ip  (2  nt  —  kn)Mb(t)dt 


rp{2 nt  -  kn)  {Mb{t)  -  Mb(kn 2~n))  d t 
■  . . .  M<;> .  .  (5.17) 


For  convenience  let  us  introduce 


A n,k(Mb)  ■=  [  il>(t)i/"’k)(dt)  and  AUjk(Mb)  :=  [  'ip{t)M(hn'k\t)dt.  (5.18) 

Jo  Jo 

Lemma  5.4.  Let  ip  be  a  wavelet  supported  on  [0,1],  and  assume  that  the  binomial 
measure  pb  satisfies  (i)-(iii)  (see  (5.5)).  Then, 

C„M  =  2“/2  .  M<;> . . .  M<;>  ■  A* (ft),  (5.19) 

and  if  ip  is  admissible  (J  ip  =  0) 

C„.k,(Mb)  =  (5.20) 

Moreover,  for  both,  fi,  and  Mb  the  random  variables  An  bn  and  A  If '  if  —  I .... ,  n)  are 
mutually  independent,  and 

An,kn  =  A0io  =  C*o,o  (5-21) 

Proof 

Use  (2.8),  (5.14)  and  (5.17)  to  obtain  the  formulas.  Independence  follows  from  the 
nested  independence  (5.6).  Finally,  the  identical  distributions  of  the  multipliers  (5.7) 
imply  that  /ip(t)p(fL'kn\dt)  is  equal  in  distribution  to  C0,o  =  Jq  fi’{t)phb{dt).  <0 

As  a  particular  example  we  find  for  the  Haar  wavelet  (use  (5.13)) 


A. 


n,kr 


=  Mr 


(n+1) 


2  kn 


Mr 


(n+1) 


2  kn-\-l 


=  M0 


Mi  —  Cq  ft. 


The  deterministic  version  has  also  been  observed  in  [7]. 

It  is  obvious  that  the  dyadic  structure  present  in  both,  the  construction  of  the  bi¬ 
nomial  measure  as  well  as  in  the  wavelet  transform,  are  responsible  for  the  simplicity 
of  the  computation  above.  It  is,  however,  standard  by  now  to  extend  the  procedure  to 
more  general  multi-nomial  cascades  such  as  Mc  introduced  in  Section  5.6  (see  [6,  91]). 
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5.4  Multifractal  Analysis  of  the  Binomial  Cascade 

Almost  all  of  the  ‘classical’  scaling  exponents  h^\  ctjA*  and  can  be  dealt  with 
explicitly  for  the  cascades.  Lemma  5.4  provides  uqf'1  readily  for  both,  Aib  and  jxb. 
Through  (5.2)  one  easily  finds  and  eventually  for  Ah;  but  increments  as  they 
appear  in  oyn')  and  h(fl>  are  not  well  defined  for  the  distribution  pb. 

For  the  sake  of  simplicity  and  to  retrace  history  let  us  start  with  the  well  known 
multifractal  analysis  of  Ah  based  on  increments.  We  will  return  to  the  wavelet  analysis 
in  Section  6.2.  As  a  starter  we  establish  the  good  news  that  for  increasing  processes 
the  partition  functions  obtained  using  or  hfj^  coincide. 

Lemma  5.5.  Assume  that  M.  is  almost  surely  increasing.  Then,  for  all  q 


Ta{q )  =  Th(q)  (5.22) 

and  for  almost  all  paths 

Ta(q,u)  =  rh(q,Lu).  (5.23) 

Proof 

To  understand  the  proof  note  that  the  task  of  estimating  partition  sums  against  each 
other  is  two-fold:  estimate  the  terms  of  the  sums  themselves,  and  ensure  that  each  term 
of  one  sum  has  a  counterpart  in  the  other  sum. 

Consider  an  increasing  path  A t(t).  For  notational  simplicity  we  set 

A[n)  =  M  ((k  +  l)2~n)  -  M(k2~n)  and  =  M  ((k  +  2)2~n)  -  M((k  -  l)2~n). 

Recall  that  we  have  =  —  (1/n)  log2  Aj^,  while  a =  —  (1  fn)  log2  Aj^. 

The  easy  part  of  the  proof  is  to  note  that  A["')  <  A[n'1  by  monotony  of  AL  Since 
there  is  a  one-to-one  correspondence  of  terms  we  find,  by  summing  over  all  k  that 
(q)  <  S^\q)  for  q  <  0  and  S^\q)  <  Sjf  -1  (q)  for  q  >  0.  Taking  limits  the  estimates 
translate  immediately  into  Th(q )  >  Ta(q)  for  q  <  0  and  ra(q)  >  Tb(q)  for  q  >  0.  By 
taking  first  expectation  and  then  limits  we  obtain  inequalities  between  the  envelopes 

n  «)■ 

To  get  inequalities  in  the  other  direction  (and  thus  establish  equality)  we  claim  first 
that  Sq1)  (q)  <  S^+2\q)  for  q  <  0.  The  difference  of  the  indices  n  and  n  +  2  will  not 
matter  in  the  limit.  For  each  ijf1^  there  is  m{k )  such  that 

[(m(k)  -  l)2~n_2,  (m(k)  +  2)2"n"2]  c  I(kn\  . 

Obviously.  <  A(fl>  again  by  monotony.  Since  rn(k)  ^  m(k')  when  k  ^  k'  we  get 

for  q  <  0  that  S^\q)  <  <  S^+2\q),  as  claimed. 

The  last  inequality  is  strict  since  there  are  only  2n  indices  of  the  form  m(k),  while 
Sf+2)  (q)  runs  over  the  entire  range  m  —  0, . . . ,  2n+2  —  1.  For  this  reason,  the  above 
argument  does  not  help  for  q  >  0.  To  obtain  an  estimate  for  q  >  0  we  need  to  construct 
a  ‘counter-term’  A(fl>  for  each  Ay"+2'1 .  But  each  interval  of  the  form  [(m  — 1)2-”-2,  (m  + 
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2)2  n  2]  intersects  exactly  two  /^"’-s.  Let  us  index  the  one  with  the  larger  increment 
by  k(m).  We  have  then  A™+2^  <  2  A  ["A  by  monotony  of  M..  Clearly,  k(m)  =  k(m') 

is  possible  even  for  m  ^  m! .  But  each  ijf1'1  intersects  only  6  intervals  of  the  form 
[( m  —  l)2_n_2,  (m  +  2)2~n^2]  and  at  the  most  6  different  values  of  m  can  lead  to  the 
same  value  k(m).  This  yields  Sln+2\q)  <  A  6  •  2q  ■  Sa  \q )  for  q  >  0, 

completing  the  proof.  <0> 

After  this  preparation  let  us  now  compute  the  multifractal  envelope  of  a  binomial 
cascade  for  a^l\  Recall  that  M.b((kn  +  l)2-n)  —  M.b(kn)  =  /Xb(lj^).  Using  in  sequence 
(5.1)  and  (5.10),  the  nested  independence  (ii)  (see  (5.6))  and  identical  distributions 
(5.7)  we  obtain 

2n— 1 

E[sth(9)]  =  E  E[(M£’)‘''  ■  ■(<’)]  EK^-tl))"] 

kn= 0 

n 

=  E[(Vlt(l))a]  ■  E  (")E  M‘  E 

i= 0 

=  ^[(M^l))9]  •  (E[Mq]  +  E[Mf])n  .  (5.24) 

Recall  that  M.b{  1)  =  almost  surely  if  the  strict  conservation  (5.5)  holds. 

Theorem  5.6.  Assume  that  A4b  is  a  binomial  cascade  for  which  (%’),  (ii)  and  (Hi) 
hold  ( see  (5.8)).  Then,  T(q)  =  —oo  if  either  o/E[Mq],  E[M^]  or  ^[(^^(l))9]  is  oo. 
Otherwise, 

TaMb(q)  =  ThMb(q)  =  -log2E  [(M0)q  +  {M1)q].  (5.25) 

Assume  furthermore,  that  the  multipliers  Ado  and  Mi  have  at  least  some  finite  moment 
of  negative  order.  Then,  when  using  increment  exponents  af1'1 ,  the  following  holds  with 
probability  one  for  all  a  such  that  T*(a )  >  0: 

dim  (K[a])  =  dim(U[al)  =  1(a)  =  f(a)  =  r*(a)  =  F(a)  =  T*(a).  (5.26) 

Assume  the  stronger  condition,  that  the  multipliers  M0  and  Mi  are  actually  bounded 
away  from  zero.  Then,  when  using  increment  exponents  tif  '1 ,  for  any  a  with  T*(a)  >  0 
the  same  equality  (5.26)  holds  almost  surely. 

The  theorem  and  (5.26)  assert  -  in  the  first  case-  that  almost  all  paths  possess  the 
same  entire  spectra,  while  in  the  second  case  the  spectra  can  be  guaranteed  to  be  equal 
only  in  a  fixed,  countable  set  of  a-s.  When  interested  only  in  a  finite  set  of  as,  however, 
this  result  is  obviously  strong  enough. 

Remark  5.7.  (Wavelet  analysis  and  fractal  support)  In  essence  the  same  formu¬ 
las  (5.25)  and  (5.26)  determine  the  spectra  emerging  from  an  analysis  of  Mb  through 
certain  analyzing  wavelet  exponents  (see  Section  6.2). 

There  is  a  natural  extension  of  this  theorem  to  multi-nomial  cascades;  however,  for 
multi-nomial  cascades  with  fractal  support  the  scaling  exponents  and  have  to 
be  adapted  [91]  to  yield  meaningful  partition  functions  (see  Section  5.6  and  (5.38)).  & 
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Remark  5.8.  (Degenerate  cascades)  Recall  that  binomial  cascades  with  conserva¬ 
tion  in  the  mean  (see  i’)  degenerate,  i.e.,  E[.Mf,(l)]  =  0,  exactly  when  T'a{  1)  <  0  (see 
[56]).  Since 

T\  1)  =  — IE[Mo  log2(M0)  +  Mi  log2(M1)]  (5.27) 

this  condition  means  that  ‘on  the  carrier’  we  multiply  mass  per  iteration  step  ‘in  the 
average’  with  a  number  smaller  than  1,  and  so  the  cascade  ‘dies  out’.  The  expression 
(5.27)  for  T'(l)  justifies  the  name  ‘information  dimension’  (compare  remark  4.6  and 
Section  5.2).  £ 

Proof 

Formula  (5.25)  is  already  established.  Corollary  4.7  shows  that  F(a)  =  T*(a).  The 
hard  part  of  the  theorem  is,  of  course,  to  find  a  formula  for  dim(iC^).  Assuming  finite 
moments  of  negative  order  for  the  multipliers  dim(iC^)  using  the  scaling  exponents 
ay"'1  was  computed  in  [8].  It  is  easy  to  verify  that  the  formula  of  [8]  for  dim(iC^) 
coincides  with  T*(a )  as  given  in  (5.25).  Earlier  results  such  as  [33]  (on  ay"'*)  and  [6] 
(on  Iipl> )  obtained  the  same  formula.  These  papers  use  the  more  restrictive  assumption 
that  the  multipliers  are  bounded  away  from  zero,  but  they  are  somewhat  easier  to  read. 
Corollary  4.1  implies  (5.26)  and  the  proof  is  complete.  <0 


5.5  Examples 

Example  5.1.  (/3eta  Binomial)  In  our  first  example  we  choose  the  multipliers  M0 
and  Mi  to  follow  a  f3  distribution,  which  has  the  density  cptp_1(l  —  t)p_1  for  t  G  [0, 1] 
and  0  elsewhere.  Hereby,  p  >  0  is  a  parameter  and  cp  is  a  normalization  constant.  Note 
that  the  conservation  of  mass  (5.5)  imposes  a  symmetrical  distribution  once  we  decide 
to  have  M0  and  Mi  —  1  —  M0  equally  distributed. 

The  /3  distribution  has  finite  moments  of  order  q  >  —p  which  can  be  expressed 
explicitly  using  the  T-function.  We  get 


/3-Binomial:  Ta(q)  = 


T(p  +  q)T(2p) 
’2Y(2p  +  q)Y(p) 


(5.28) 


and  T(q)  =  — oo  for  q  <  —p.  For  a  typical  shape  of  these  spectra  see  Figure  6. 

In  [89]  binomial  cascades  were  successfully  used  for  modelling  data  traffic  on  the 
internet.  To  provide  more  flexibility  the  distributions  of  the  multipliers  were  allowed 
to  depend  on  the  scale  n  (compare  Example  7.2).  4k 

Example  5.2.  (Uniform  Binomial)  As  a  special  case  of  the  j3e ta  binomial  we  obtain 
uniform  distributions  for  the  multipliers  when  setting  p  —  1.  The  formula  (5.28)  sim¬ 
plifies  to  Ta(q)  =  —  1  +  log2(l  +  q)  for  q  >  —  1.  Applying  the  formula  for  the  Legendre 
transform  (4.11)  yields  the  explicit  expression 
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FIGURE  6.  The  spectrum  of  a  binomial  measure  with  (3  distributed  multipliers  with  p  =  1.66. 
The  formula  is  given  by  (5.28).  Trivially,  T( 0)  =  —1.  whence  the  maximum  of  T*  is  1.  In 
addition,  every  positive  increment  process  has  T(  1)  =  0,  whence  T*  touches  the  bisector. 
Finally,  the  LRD  parameter  is  fhvar  =  (T( 2)  +  l)/2  =  0.85)  (see  (7.20)  and  (7.25)  below). 


Example  5.3.  (Log-normal  Binomial)  Another  case  of  strong  interest  is  the  use  of 
log-normal  distributions  for  the  multipliers  M0  and  M\ .  Note  that  we  have  to  replace 
(i)  by  (i’)  (see  (5.8))  in  this  case  since  log-normal  variables  can  be  arbitrarily  large,  i.e., 
larger  than  1.  We  comment  on  some  peculiarities  this  might  cause  after  computing  the 
spectrum. 

Deterministic  envelope:  As  the  multipliers  M0  and  M\  of  a  log-normal  binomial 
are  log-normal  variables  we  need  expressions  for  the  moments  of  such  variables.  Let 
G  be  a  Gaussian  with  mean  m  and  variance  a2,  and  set  M  =  eG .  Then,  IEfM9]  = 
IE[exp(gG)]  =  exp(gm  +  q2cr2/ 2). 

For  further  simplification  let  us  consider  the  case  where  Mo  and  Mi  are  equally 
distributed.  Then,  their  mean  must  be  1/2  which  translates  into  m  +  <j2/2  =  —  ln(2). 
We  find  from  (5.25)  that 


log-normal  binomial: 


T«(q) 


(q-1) 


2 ln(2) q ) 


(5.30) 


for  all  q  G  IR  such  that  E^At^l))9]  is  finite.  For  later  use  note  that  the  parabola  in 
(5.30)  has  two  zeros:  1  and  gcrit  =  21n(2)/a2. 

Degenerate  case  and  tail  of  marginals  It  is  known  also  for  more  general  cascades 
that  1E[A/1{)(1)]  =  0  if  and  only  if  T'(l)  <  0  [56].  Since  A4fe(l)  is  a  positive  random 
variable  this  occurs  if  and  only  if  Af;,(l)  =  0  almost  surely,  and  we  say  that  the 
cascade  is  degenerate. 

To  ensure  non-degenerate  cascades  for  the  remainder  we  assume  T'(l)  >  0.  Due  to 
the  parabolic  shape  of  (5.30)  this  condition  may  be  expressed  easily  as  gcrit  >  1,  i.e., 
u2  <  2  ln(2). 

Assuming  that  T'(  1)  >  0,  [56]  shows  that  T'(l)  gives  the  almost  sure  dimension 
of  the  carrier  of  the  cascade.  Since  paths  with  .Ad;,(l)  =  0  have  an  empty  carrier  we 
conclude  that  actually  M.h(l)  ^  0  almost  surely,  not  only  with  positive  probability. 
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Always  assuming  T'{  1)  >  0,  another  result  of  [56]  states  that  IE^Al^l))9]  <  oo 
exactly  for  q  <  gcrit.  Moreover,  [66]  conjectures  and  [75]  proves  the  marginal  tail  be¬ 
havior  P[WG(1)  >  x]  ~  x~qciit.  We  will  return  with  an  explanation  in  an  instant. 
For  the  moment  we  conclude  that  (5.30)  is  valid  exactly  for  q  <  qcrit,  for  all  other  q 
T(q)  =  -oo. 

Legendre  spectrum:  Since  T(g)  is  differentiable  exactly  for  q  <  qcr it  we  may  obtain 
its  Legendre  transform  implicitly  from  (4.11)  for  a  =  T\q )  with  q  <  qcrit,  i.e.,  for  all 
a  >  acrit  =  T'(gcrit)  =  <x2/(21n(2))  —  1.  Eliminating  q  from  (4.11)  yields  the  explicit 
form 

r»  =  1  -  (“  -  1  -  («  >  ««)■  (5.3i) 

For  a  <  acrit  the  Legendre  transform  yields  T*(a)  =  a  ■  gcrit.  Thus,  at  acrit  the  spectrum 
T*  crosses  over  from  the  parabola.  (5.31)  to  its  tangent  through  the  origin  with  slope 
gcrit  (the  other  tangent  through  the  origin  is  the  bisector). 

Since  gcrit  >  1  we  have  acrit  <  0  and  the  spectrum  possesses  negative  a  with  finite, 
though  negative  T*(a)  (compare  with  the  ‘uniform’  binomial  cascade  above  which  has 
no  such  a).  We  return  to  this  fact  in  an  instant.  The  two  zeros  of  T*  are  both  positive 
for  any  choice  of  a,  since  T*  is  — apart  from  the  linear  part  for  a  <  acrit-  a  parabola 
which  touches  the  bisector  and  has  a  maximum  value  of  1. 

Rescaled  histograms:  The  log-normal  framework  allows  also  to  calculate  F (a) 
explicitly.  We  use  the  opportunity  to  demonstrate  which  rescaling  properties  of  the 
marginal  distributions  of  the  increment  processes  of  A4&  are  captured  in  the  multifractal 
spectra. 

First,  note  that  —  ln(2)  •  aj.n)  =  (l/n)ln (nb(I^))  (see  (2.13)  and  (5.10))  is  the  sum 
of  n  independent  Gaussian  random  variables  (1/n)  ln(M of  mean  m/n  and  variance 
o2jn2,  plus  an  independent  additive  term  distributed  as  (1/n)  ln(A4f,(l)).  To  keep 
computations  simple  we  neglect  this  last  term;  indeed,  it  will  converge  in  distribution 
to  0  since  ln(A46(l))  has  a  Weibull  tail  and,  thus,  will  not  affect  the  distribution  of 
in  the  limit.  In  summary,  this  amounts  to  approximating  —  ln(2)  •  a ^  by  a  Gaussian 
of  mean  m  and  variance  a2/n.  The  mean  value  theorem  of  integration  gives 


-Pfi[|c4n)  -  a  |  C  e\  ~ 


fln(2)(— a+e) 


sj2'Ka2/n  </ln(2)(— a— e) 

1 


exp 


(x  —  mfi 


dx 


y/'2ir<j2/n 


In (2)  •  2e  •  exp 


2  o2  fin 
(-  ln(2)ga|W  -  m)2 
2  o2  jn 


with  xajn  G  [a  —  e,a  +  e]  for  all  n.  To  compute  F(a)  we  first  need  to  sum  over  k  (see 
(3.22));  however,  the  expression  is  identical  for  all  k.  Next,  keeping  only  the  exponential 
term  in  n  and  substituting  m  =  —a2/2  —  ln(2)  we  find  for  large  n 


1 

n 


(5.32) 


Taking  hnally  the  limit  of  e  — »•  0  we  see  that  F(a)  equals  the  right  hand  side  of  (5.31) 
for  all  a,  and  thus  T*(a)  =  F(a)  for  all  a  =  T'(q)  as  stated  in  corollary  4.7.  On  the 
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other  hand,  T*(a)  >  F(a)  for  a  <  acrit-  The  above  computation  shows  impressively  how 
well  adapted  a  multiplicative  iteration  with  log-normal  multipliers  is  to  the  multifractal 
analysis  (or  vice  versa):  F  extracts,  basically,  the  exponent  of  the  Gaussian  kernel. 

In  the  range  of  a  with  F(a)  >  0  — where  equality  in  the  multifractal  formalism  is 
known  to  hold  for  Mb  (see  (5.26)) —  these  features  can  be  measured  or  estimated  from 
a  ‘typical’  realization  via  the  re-normalized  histogram ,  i.e.,  the  grain  based  multifractal 
spectrum  /(a).  This  is  a  property  which  could  be  labelled  with  the  term  ergodicity. 
Note,  however,  that  classical  ergodic  theory  deals  with  observations  along  an  orbit 
of  increasing  length  while  /(a)  is  in  terms  of  a  sequence  of  orbits.  Also  recall  that 
estimating  the  negative  F(a)  requires  using  a  large  array  of  realizations  of  the  process 
(compare  remark  3.17). 

Partition  function: 

Since  almost  all  paths  of  cascades  are  increasing,  the  scaling  exponents  cq”'1  must 
eventually  become  positive  for  any  such  path  and  ra(q)  must  be  increasing.  The  de¬ 
terministic  envelope  Ta(q),  on  the  other  hand,  has  a  parabolic  shape  where  it  is  finite. 
So,  r(q)  >  T(q)  for  q  >  (1  +  gcrit)/2,  i.e.,  for  q  beyond  the  maximum  of  T(q). 

Log-normal  marginals:  The  log-normal  binomial  enjoys  the  advantage  of  hav¬ 
ing  exactly  log-normal  marginals  Hb(lj^)  since  the  product  of  independent  log-normal 
variables  is  again  a  log-normal  variable. 

Finally,  we  return  to  the  various  issues  we  have  left  unresolved  above. 

Negative  (virtual)  exponents:  The  most  puzzling  consequence  of  the  absence  of 
strict  conservation  of  mass  resides  in  the  presence  of  negative  a  with  finite  T*(a).  The 
apparent  paradox  of  such  a  was  first  mentioned  and  studied  in  [66] :  negative  exponents 
seem  to  indicate  that  increments  over  smaller  child-intervals  could  be  larger  than  the 
increment  over  the  parent-interval  which  is  impossible  for  increasing  paths. 

The  paradox  was  explained  away  in  [66]  with  the  observation  that  virtual  exponents 
should  be  interpreted  not  as  a  property  of  the  limiting  cascade  and  the  actual  process, 
but  rather  as  a  property  of  the  iterative  cascade  construction  as  the  construction  pro¬ 
gresses,  what  we  expect  an  increment  to  be  over  an  interval  might  increase  (recall  that 
(5.11)  is  not  the  actual  increment  but  its  expectation). 

Here,  we  further  demystify  virtual  exponents  in  providing  an  interpretation  in  terms 
of  the  actual  process,  not  only  its  construction.  As  we  have  seen  above  and  as  follows 
from  corollary  4.7  with  more  generality  we  have  F(a)  =  T*(a)  for  a  =  T\q).  For  such 
a  F(a)  is  thus  finite,  and  IE[Adn)(a,  e)]  >  0.  If  a  is  in  addition  virtual,  then  F(a)  <  0 
and  the  expected  number  of  cxfi'1  ~  a  decays  rather  rapidly  as  resolution  n  increases 
but  never  reaches  zero.  Thus,  even  for  a  large  n  an  increment  larger  than  1  over  an 
interval  of  length  2~n  can  be  observed,  though  extremely  rarely,  leading  to  a  negative 
.  This  does  not  contradict  the  fact  that  for  any  fixed  path  all  will  eventually 
(for  large  n)  be  positive. 

* 


Remark  5.9.  (Deterministic  envelope:  shape  and  interpretation) 

There  are  several  distinct  parts  of  the  ‘Legendre  spectrum’  T*(a). 

Looking  first  at  T* ,  the  positive  part  (where  T*(o)  >  0)  relates  most  directly  to  local 
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properties  of  individual  paths.  For  cascades,  it  provides  dim(£'^),  the  dimension  of 
sets  with  equal  local  behavior  as  measured  through  s^\  as  well  as  /(a)  and  r*(a);  in 
general  it  provides  at  least  an  upper  bound.  This  part  of  T*(a )  lies,  obviously,  between 
its  two  zeros  which  in  turn  correspond  to  the  slopes  of  the  two  tangents  at  T(q)  through 
the  origin  (recall  that  these  functions  are  concave  and  apply  the  Legendre  formalism 
(4.7)). 

The  remaining  T*(a),  i.e.,  the  negative  ones,  can  not  be  interpreted  as  actual  dimen¬ 
sions.  However,  through  the  LDP  of  (4.11)  and  (3.23)  they  still  provide  information  on 
particularly  rare  degrees  a  of  local  scaling.  Sometimes  they  are  referred  to  as  ”  negative 
dimensions” . 

Second,  let  us  discuss  now  the  values  of  a  which  may  appear  in  the  spectrum.  The 
not  immediately  intuitive  ones  are  the  ’’latent”  exponents  (a  with  T*(a)  <  0)  and  the 
”  virtual”  ones  (a  <  0).  For  the  scaling  exponents  ay”"1  and  of  increasing  processes 
all  virtual  exponents  are  necessarily  latent.  This  follows  from  T(l)  =  0  which  implies 
that  T*(a)  <  a  (compare  remark  4.6);  indeed,  no  increasing  path  can  show  a  negative 
increment-exponent  in  the  limit  and  virtual  exponents  can  exist  here  only  as  ‘rare 
events’  (compare  remark  3.17).  The  name  ‘virtual’  has  to  be  understood  in  this  context. 
For  wavelet-based  exponents  w[n'>  and/or  non-increasing  processes,  however,  negative 
a  may  occur  pathwise  and  very  naturally  (compare  Section  6.1). 

* 


5.6  Multiplicative  Schemes  beyond  Dyadic  Structure 

The  dyadic  structure  inherent  to  the  binomial  cascades  can  at  times  be  too  restrictive. 
The  broader  class  of  ‘self-similar  cascades’  offers  more  flexibility  going  beyond  dyadic 
and  even  towards  random  subdivision  of  intervals.  Despite  this  added  freedom  one  can 
think  of  cascades  in  principle  as  various  types  of  ‘muscles’  hanging  off  an  always  present 
tree-skeleton.  Indeed,  the  tree  at  the  core  of  all  cascades  provides  common  structure 
and  a  central  tool  in  the  multifractal  analysis,  enabling  an  approach  similar  to  the 
binomial  case. 

It  should  not  go  without  mentioning,  though,  that  various  technical  difficulties  arise 
when  dealing  with  this  more  general  framework,  such  as  with  the  moments  of  negative 
order  for  cascades  with  fractal  support.  With  regard  to  this  particular  problem  we 
describe  here  a  method  which  relies  on  modifying  the  scaling  exponents  as  in  [91]. 
For  wavelet  exponents,  the  way  to  go  is  to  modify  the  partition  function  by  employing 
the  lines  of  maxima,  pioneered  in  [78].  It  is  notable  that  both  approaches  can  be 
considered  as  a  multifractal  analysis  with  new  scaling  exponents,  related  to  and 
Wjf*  respectively;  consequently,  the  multifractal  formalism  of  section  4  still  holds  for 
these  modified  methods. 

Example  5.4.  (Multi-nomial  cascades)  The  first  most  obvious  generalization  of 
the  binomial  measure  is  to  split  intervals  into  c  subintervals.  Making  the  obvious  ad¬ 
justments  such  as  organizing  the  multipliers  now  in  a  c-ary  tree  and  adapting  (i)-(iii) 
(see  (5.5))  this  leads  to  a  measure  nc. 
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Towards  understanding  the  multifractal  properties  of  /ic  choose  m  such  that  2n  ~  cm. 
Then,  a  typical  interval  of  length  2~n  will  have  a  /ic-rnass  approximately  equal  to 
M(S0)  ■  ■  ■  ■  M where  k,  indicates  now  position  on  a  c-ary  tree.  As  in  (5.25)  we 

find  IEa[5'^nj/(c(g)]  ~  (IE[(M0)9  +  . . .  +  (Mc_i)q])m.  Since  m/n  ~  log(2)/ log(c)  this 
leads  to 

TaMM  =  -  logcIE[(Mo)*  +  . . .  +  (M^Y]  (5.33) 

This  formula  is  indeed  correct  for  all  q  G  IR  and  equality  holds  in  the  multifractal 
formalism  provided  all  the  multipliers  are  zero  with  probability  zero. 

Cascades  with  fractal  support  and  scaling  exponents:  A  true  advantage  of 
splitting  intervals  into  more  than  two  subintervals  resides  in  the  possibility  to  leave  some 
of  these  subintervals  without  mass:  split  intervals  into  d  >  c  subintervals  and  distribute 
mass  in  the  usual  multiplicative  manner  but  only  among  selected  c  subintervals.  This 
construction  can  be  viewed  as  an  extreme  case  of  a  c'-nomial  cascade  where  now  certain 
multipliers  are  set  to  zero  almost  surely.  This  leads  to  conjecture  that 

TmM )  =  -  logc/  E[(M0)«  +  . . .  +  (Mc-i)q]  (5.34) 

Due  to  the  construction  the  limiting  measure  has  a  fractal  support  of  dimension 
— T(0)  =  log(c)/log(c/)  <  1.  This  might  be  of  practical  advantage:  when  modelling 
rain  data,  e.g.,  the  actual  events  of  interest  occupy  a  very  small  set  on  the  time  axis. 

It  may  be  ’shocking’  news  at  first  sight  that  equality  in  the  multifractal  formalism 
in  the  sense  of  theorem  5.6  holds  no  longer  for  cascades  with  fractal  support  (see  [91]). 
What  breaks  down  is  the  partition  function:  Unless  c  is  a  power  of  2,  (5.34)  holds  for 
exponents  and  a ^  only  for  q  >  0  and 

Ta(q)  =  Ta(q )  =  -oo  ( q  <  0)  (5.35) 


due  to  boundary  effects  which  produce  exponents  that  are  exceptionally  large,  yet  finite. 
As  is  demonstrated  in  [91]  (see  also  [84]),  a  simple  remedy  for  increasing  processes  AA 
is  to  exclude  the  dyadic  intervals  fj”'1  over  which  M.  is  constant  and  to  set 


h ^  if  a ^  y  oo 
oo  otherwise. 


(5.36) 


(n) 


is  the  increment  of  M.  over  the  interval  U  I ^  U  lj^1  if  M. 


In  other  words,  2-nAfc 

has  a  strictly  positive  increase  over  the  middle  interval  ij!''1 :  if  A4  is  flat  over  the  middle 
interval  then  2_"Afc  )  =  0  and  there  is  no  contribution  of  I ^  to  the  partition  sum. 

Let  us  check  that  this  procedure  is  an  extension  of  the  usual  analysis.  First,  since 
limits  of  /q"'1  and  (n  — >  oo)  are  equal  we  have  E\f  =  etc.  Second,  if  M-  is 
strictly  increasing  then  h ^  =  A)"-  even  for  all  n  and  /c,  and  T\  =  tv  Finally,  for  general 
cascades  and  positive  q  (which  correspond  to  small  a  of  the  spectrum)  the  partition 
function  has  not  changed  either: 

Lemma  5.10.  [91]  Assume  that  M.  is  increasing.  Then, 


T\ (q)  =  rh(q)  =  ra(q)  (q  >  0). 


(5.37) 
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To  obtain  the  multifractal  properties  of  Aic  take  expectations  in  the  proofs  of  [91].  It 
follows  that  Tx,Mc(q)  is  the  unique  solution  of  (5.34)  for  all  q.  From  [33,  6,  8]  it  follows 
then  that  the  dimension  based  spectra  are  equal  to  T£  where  it  is  positive.  In  short, 
equality  holds  in  the  multifractal  formalism  in  the  following  sense:  For  multi-nomial 
cascades 

dim(Aja])  =  dim(Ajal)  =  fx,Mc(a)  =  Tx,Mc(a )  =  T\Mc(a )  (5-38) 

wherever  the  right  hand  side  is  positive. 

As  is  shown  in  [84,  91]  using  the  coarse  singularity  exponent  A[nl  provides  a  reliable 
estimation  of  the  multifractal  structure  of  distributions  with  fractal  support.  Using 
or  aj,n\  on  the  other  hand,  leads  to  numerical  instability  when  estimating  r(q)  for 
negative  q.  4 

A  more  radical  way  of  leaving  the  binary  geometric  structure  behind  is  by  random¬ 
izing  it.  To  this  end  we  consider  first  the  inverse  of  binomials.  Inverse  measures  can  be 
defined  in  all  generality  (see  [72,  94])  allowing  for  further  flexibility. 

Example  5.5.  (Random  geometry:  Inverse  measures)  Let  Ai\  be  the  inverse 
function  of  a  binomial  distribution  function  Aib,  in  other  words,  Aib(Aib(t))  =  t. 

Then,  Ail  is  almost  surely  increasing  and  continuous.  So,  it  defines  a  random  measure 
/i|  which  can  be  described  as  follows:  for  all  n,  the  unit  interval  is  divided  into  2” 
subintervals  J^l)  of  random  lengths  and  with  /i\  mass  1/2”  each. 

To  motivate  a  formula  for  T(q )  in  this  case  let  us  note  that  2”7  ■  IEfS'^^g)]  should 
tend  to  0  for  7  >  T(q),  resp.  to  oo  for  7  <  T(q).  Noting  that  2”7  =  |/[n)|~7,  we  may 
equivalently  study  E[J0 /u^(J^m))9|./^u)|~7]  for  hirer  and  finer  partitions  of  [0,1]  into 
intervals  J^n> .  Therefore,  we  iteratively  split  the  largest  interval  of  such  a  partition, 
say  jiu\  into  its  subintervals  J^+V>  and  which  have  lengths  M scaleuo2V  •  \JvU)\ 

and  •  |  Jlu\-  The  contribution  of  J ^  to  the  partition  sum  was  Ab(  Jv^)q\  Jv^  |-7; 

the  iteration  step  will  replace  it  by  Afc('^2v+1^)9l'4^+i'l_7  +  A^^+i^l^+i^l-7  which 
can  be  written  as 

dOtTlU’n  (2-a(M<;+11)-^  +  2-«(A4"i11)-/  .  (5.39) 

This  leads  us  to  pose  that  TM\  (q)  is  the  unique  solution  of 

E[(1/2)«(M0)-t(9)  +  (1  /2)q(M1)~T^]  =  1.  (5.40) 

Indeed,  during  iteration  IE[^/) /r^(  J^)9|  J^|-7]  will  decay,  resp.  increase  (on  the  av¬ 
erage)  according  to  whether  7  >  T(q)  or  7  <  T(q),  as  required.  4 

Example  5.6.  (Statistically  self-similar  measures)  Combining  above  examples 
leads  to  defining  a  random  binomial  Aic  by  splitting  intervals  iteratively  into  c 
subintervals  Jc(/+1) , . . . ,  J^Xc-i  with  lenSth  I  4fc+*1)  I  =  LS+i  ’  1 4" ’  I  and  mass  )  = 

M-ck+P  hi  the  most  simple  case,  one  will  require  conservation  of  mass,  i.e., 
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Mj+1)  +  . . .  +  =  1,  but  also  L^+1)  +  . . .  +  Z/^fc+c- 1  =  1  which  guarantees  that 

/jc  lives  everywhere.  Assuming  the  conditions  (ii)  and  (iii)  (see  (5.7))  hold  for  both,  the 
length- multipliers  L^l>  as  well  as  the  mass-multipliers  we  hnd  that  Tj^c(q)  is  the 
unique  solution  of 

E[(M0)9(L0)-t(9)  +  . . .  +  (Mc_1)q(Lc_i)-T^]  =  1.  (5.41) 

A  rigorous  proof  of  (5.41)  is  obtained  by  ‘taking  expectations’  in  the  proof  of  [91,  Prop 
14],  Doing  so  shows,  moreover,  that  T(q)  assumes  a  limit  in  these  examples.  4k 

In  the  broader  context  of  Example  5.6  it  is  easy  to  recognize  that  all  the  examples 
above  are  merely  special  cases  of  (5.41).  This  confirms,  in  particular,  (5.40). 

Example  5.7.  (Cascades  in  the  wavelet  domain)  As  the  concluding  example  we 
mention  that,  with  regard  to  (5.14),  one  may  choose  to  model  directly  the  wavelet  coef¬ 
ficients  of  a  process  in  a  multiplicative  fashion  in  order  to  obtain  a  desired  multifractal 
structure.  First  steps  in  this  direction  have  been  taken  in  [1],  4k 

A  final  remark  is  in  order. 

Remark  5.11.  (Stationary  increments:)  Due  to  property  (iii)  the  binomial  cascade 
has  increments  A4  (( k  +  l)2~n )  —  Ai(k2~n)  which  have  identical  distributions  provided 
Mq  and  M\  are  equally  distributed.  This  first  order  stationarity  of  dyadic  increments  is 
sometimes  sufficient  to  obtain  interesting  results  (see  lemma  8.8).  On  the  other  hand, 
these  increments  are  clearly  not  second  order  stationary. 

However,  to  obtain  true  stationarity  of  increments  one  has  to  leave  rigid  multi¬ 
nomial  subdivisions  towards  stationary  ones,  e.g.,  produced  by  Poisson  arrivals.  For 
some  appealing  examples  see  [14,  70,  74, 105].  With  these  models  one  trades,  potentially, 
appealing  statistical  properties  against  more  involved  analysis.  & 


6  Wavelet  based  Analysis 

In  this  section  we  start  with  a  note  on  the  appealing  connection  between  wavelet- 
based  multifractal  analysis  and  classification  of  functions  using  regularity  spaces  such 
as  Besov  spaces.  Then,  we  revisit  the  binomial  cascades  3with  wavelets  and  demonstrate 
a  simple  relation  between  the  multifractal  properties  of  a  process  and  its  (distributional) 
derivative. 


6.1  Wavelet  Coefficients  and  Besov  Spaces 

Besov  spaces  are  function  spaces  tailored  for  studying  global  regularity.  They  have 
become  of  fashion  recently  since  an  elegant  description  in  terms  of  wavelet  coefficients 
(see  [76])  has  become  available.  In  [76]  it  is  shown  that  the  norm  of  the  Besov  space 
B°(LU)  of  a  process  with  wavelet  coefficients  Ch\.  is  equivalent,  in  the  notation  of  (2.7), 


6  Wavelet  based  Analysis  51 


to 


Roughly  speaking,  this  norm  measures  the  smoothness  of  order  s  in  Lu,  where  v  is  an 
additional  parameter  for  making  finer  distinctions  in  smoothness. 

A  connection  between  multifractal  theory  and  Besov  spaces  has  been  established  in 
[30].  Here,  we  note  that  multifractal  analysis  (using  convenient  orthogonal  wavelets) 
can  be  viewed  as  determining  in  which  Besov  spaces  the  analyzed  process  lies  [89]. 
Indeed,  by  the  very  definition  (3.10)  of  the  partition  function  with  (see 

(2.11))  the  Bl(Lu)  norm  of  a  path  of  the  process  is  finite  if  su  <  rw(u)  +  1  and  infinite 
if  su  >  tw(u)  +  1  (see  [53]  for  similar  results).  In  other  words,  pathwise 

sup{s  :  Y  G  BSV(LU)}  =  (6.2) 

Strictly  speaking  s  must  be  smaller  than  the  regularity  r  of  the  wavelet  ip  for  (6.1) 
and  (6.2)  to  hold.  In  other  words,  at  least  r  vanishing  moments  as  well  as  r  continuous 
derivatives  of  ip  are  required.  Given  this,  Besov  norms  do  not  depend  on  the  choice  of 
the  wavelet  basis  up  to  bounded  factors.  Since  the  multifractal  analysis  using  wavelets 
determines  the  Besov  spaces  that  contain  the  signal,  we  conclude  that  rw(u)  will  not 
depend  on  the  choice  of  the  wavelet,  provided  the  above  regularity  conditions  are  met. 

6.2  The  Binomial  Revisited  with  Wavelets 

The  binomial  measure  pb  is  not  a  function  or  process  in  the  usual  sense,  but  a  random 
distribution.  Thus,  increment  exponents  such  as  /q"'1  and  oqn')  are  not  well  defined  for  pb 
and  we  have  to  rely  on  wavelet  exponents  u/fl>  to  reveal  the  local  scaling  structure.  As 
is  often  the  case,  the  deterministic  envelope  is  the  most  simple  spectrum  to  compute. 

Taking  into  account  the  normalization  factors  in  (2.11)  when  using  lemma  5.4,  the 

calculation  of  (5.24)  carries  over  to  give 

ElSilto]  =  2”’E[|C0,„|«]  •  (EnM  +  En[M?])“ . 

Similar  manipulations  work  for  S^'M  .  In  summary: 

Lemma  6.1.  Provided  the  wavelet  coefficients  !E[|C'o)o(/-t&)|9]  resp.  IE[|C,o,o(AI{,)|'?]  and 
the  moments  IE [Mq],  lE[Mf]  and  IE[(.Mb(l))9]  are  all  finite  we  have 

Tw,i,b(q)+q  =  TwMb(q)  =  TaMb(q),  (6.3) 

Tw,»bia^1)  =  Tw,mM)  =Ta,Mb(a)-  (6-4) 

Imposing  additional  assumptions  on  the  distributions  of  the  multipliers  we  may  also 
control  v/jf  (pf)  themselves  and  not  only  their  moments.  To  this  end,  we  should  be  able 
to  guarantee  that  the  (Haar)  wavelet  coefficients  don’t  decay  too  fast  (compare  (2.10)), 
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i.e.  the  random  prefactor  RHS  in  (5.13)  and  (5.14)  doesn’t  become  too  small.  Indeed,  it 
is  sufficient  to  assume  that  there  is  some  £  >  0  such  that  | Co,o  I  >  e  almost  surely,  (recall 
that  Co,o  =  2M0  — 1  for  the  Haar  wavelet.)  Then  for  all  t,  (1  /n)  log (f  if n^’kn\dt))  — >  0, 
and  with  (5.14) 

wM6(t)  =  liminf  --log2  (2n/2|C„ifcfl|)  =  aMh(t)  ~  1,  •  (6.5) 

n—>oc  Tl 


This  relation  agrees  with  (6.4)  and  appeals  to  one’s  intuition:  the  singularity  exponents 
of  a  ‘sufficiently  nice’  process  and  its  (distributional)  derivative  should  differ  by  1,  as 
is  demonstrated  here  with  binomial  cascades.  We  should  point  out,  though,  that  this 
simple  rule-of-thumb  may  fail  miserably  in  the  presence  of  ‘chirps’,  even  for  increasing 
processes.  More  precise  information  on  how  scaling  properties  change  locally  — such  as 
provided  by  the  ‘two-microlocalization’  analysis  [51,  50] —  is  required  to  fully  under¬ 
stand  how  taking  derivatives  affects  the  local  Holder  regularity  (see  also  Section  6.3). 

Corollary  6.2.  Assume  that  fib  is  a  binomial  measure  satisfying  (i)-(iii)  (see  (5.5)). 
Assume  furthermore,  that  the  random  variables  An^k(pb)  resp.  Ank(M.b )  are  uniformly 
bounded  away  from  0.  Then,  equality  holds  in  the  multifractal  formalism  for  the  wavelet 
based  spectra  of  pb,  resp.  A4b,  meaning  that 


TatMb(a  +  !)’ 


(6.6) 


respectively 

=  /■».*.(<>)  =  <MbM  =  W»).  (6-7) 

For  deterministic  binomials  twm  has  been  calculated  in  [7].  Note  that  for  the  Haar 
wavelet  An,k  =  M2(£+1)  -  =  2 M%+1)  -  1. 

Requiring  that  Ank  should  be  bounded  away  from  zero  in  order  to  ensure  (6.5), 
though  satisfied  in  some  simple  cases,  seems  unrealistically  restrictive  to  be  of  prac¬ 
tical  use.  A  few  remarks  are  in  order,  then.  First,  this  condition  can  be  weakened  to 
allow  arbitrarily  small  values  of  Anj.  as  long  as  all  its  negative  moments  exist.  This 
can  be  shown  by  an  argument  using  the  Borel-Cantclli  lemma.  Furthermore,  if  one 
drops  (iii)  and  allows  the  distributions  of  the  multipliers  to  depend  on  scale  (compare 
Example  7.2),  then  An  k  has  to  be  bounded  away  from  zero  only  for  large  n.  In  applica¬ 
tions  such  as  network  traffic  modelling  one  finds  indeed  that  on  fine  scales  An  k  is  best 
modelled  by  discrete  distributions  on  [—1,1]  with  large  variance  and  no  mass  around 
0. 

Another  way  out  is  to  avoid  small  wavelet  coefficients  at  all  in  a  multifractal  analysis. 
More  precisely,  one  would  follow  [7,  53]  and  replace  Cn^n  in  the  definition  of  v/f  '1  (2.11) 
by  the  maximum  over  certain  wavelet  coefficients  ‘close’  to  t.  Since  this  procedure 
amounts  to  nothing  more  than  introducing  a  yet  better  adapted  scaling  exponent  w , 
the  multifractal  formalism  still  holds  in  the  sense  of  Section  4.  Towards  establishing 
equality  between  spectra  [53]  gives  conditions  under  which  the  spectrum  (a)  based 
on  this  new  exponent  w  agrees  with  the  ‘classical  Holder  exponent  spectrum’  dim(S^) 
based  on  h^‘\A4b). 
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6.3  Multifractal  Properties  of  the  Derivative 

Corollary  6.2  establishes  for  the  binomial  cascade  what  intuitively  appears  to  be  of 
general  truth:  taking  derivatives  should  decrease  the  local  Holder  regularity  by  1.  This 
section  is  devoted  to  a  few  words  of  warning  in  this  context.  Let  us  start  with  a  simple 
example  where  the  given  intuitive  rule  of  thumb  fails  utterly  and  continue  then  to 
discuss  a  practical  but  powerless  approach  to  studying  the  multifractal  properties  of 
the  derivative. 

As  indicated  in  the  last  section,  the  above  simplistic  rule  of  thumb  fails  in  the  presence 
of  chirps  [51].  For  an  example  consider  the  function  M(t)  =  5 t  +  f3sin(l/f).  Since  its 
derivative  M. '  is  strictly  positive,  this  is  a  strictly  increasing  function.  Furthermore,  the 
local  Holder  regularity  of  M.  at  t  —  0  is  of  order  3,  i.e.,  Hm( 0)  =  3,  since  | M.(t)  —  5t\  < 
|f|3  is  the  best  polynomial  approximation  in  the  sense  of  (2.1).  Its  derivative,  however, 
shows  HMf{ 0)  =  1  since  it  can  be  approximated  as  5|  <  \t\.  The  local  regularity 

dropped,  thus,  by  2.  This  is  due  to  the  strong  oscillations,  also  called  a  chirp,  present 
at  t  —  0.  To  fully  understand  and  assess  how  the  differential  operator  affects  the  local 
Holder  regularity,  one  has  to  get  acquainted  with  the  ‘two-microlocalizatioiT  analysis 
[51,  50], 

Now,  to  explore  the  role  of  wavelets  in  this  matter  we  place  ourselves  in  the  framework 
where  one  cares  less  for  a  wavelet  representation  but  where  one’s  interests  he  solely 
in  an  analysis  of  the  oscillatory  behavior.  Therefore,  we  will  employ  analyzing  mother 
wavelets  0  which  are  typically  not  orthogonal,  such  as  the  derivatives  of  the  Gaussian 
kernel  exp  (— f2/ 2)  which  were  used  to  produce  Figure  7. 

To  connect  wavelet  coefficients  of  a  process  and  its  derivative  we  use  integration  by 
parts.  For  a  continuous  measure  //  on  [0, 1]  with  distribution  function  A4(t)  =  /x([0,  t)) 
and  a  continuously  differentiable  function  g  this  reads  as 
r  2n  — 1 

g(t)n(dt)  =  hm  s(fe2-”)  (A^((fe  +  1)2-")  -  A^(fe2-”)) 

J  rwoo 

2n  — 1 

=  lim  V  M(k2~")  ( g((k  -  1)2-”)  -  9(KT")) 

n—>oc  z ' 
k= 0 

+M(l)g(l-2-n)-M(0)g(-2~n) 

=  M(l)g(l)  -  M(0)g(0)  -  J  M{t)g\t)dt  (6.8) 

where  we  alluded  to  (5.3)  and  regrouped  terms.  If  M.  has  a  derivative,  then  /i(dt)  be¬ 
comes  M.'{t)dt  and  (6.8)  reduces  to  the  usual  rule  of  integration  by  parts.  Now,  setting 
git)  =  2n/2'0(2 nt  —  k)  for  a  smooth  analyzing  wavelet  ifj  we  have  g'it)  =  23n/2'0/(2nf  —  k) 
and  using  that  M.{1)  —  1  and  Ad(0)  =  0,  we  get 

Cn,k(fp,fk)  =  2^(2”  -  k)  -  2n  ■  Cn^',M).  (6.9) 

Estimating  2”  —  kn  =  2n  —  \t2n\  ~  (1  —  t) 2n  and  assuming  exponential  decay  of  -0(f) 
at  infinity  allows  to  conclude 


=  -1  + 


(6.10) 
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FIGURE  7.  Demonstration  of  the  multifractal  behavior  of  a  binomial  measure  pb  (left)  and 
its  distribution  function  Mb  (right).  On  the  top  a  numerical  simulation,  i.e.,  (5.2)  on  the  left 
and  Mb(k2~n )  on  the  right  for  n  =  20.  In  the  middle  the  moduli  of  a  continuous  wavelet 
transform  [23]  where  the  second  derivative  of  the  Gaussian  was  taken  as  the  analyzing  wavelet 
ip(t)  for  pb,  resp.  the  third  derivative  ip'  for  Mb-  The  dark  lines  indicate  the  ‘lines  of  maxima’ 
[53,  7],  i.e.,  the  locations  where  the  modulus  of  f  ipplH.  —  s)pb(  di)  has  a  local  maximum  as  a 
function  of  s  with  j  fixed.  On  the  bottom  a  multifractal  analysis  in  three  steps.  First,  a  plot 
of  log  Sy)'1  (q)  against  n  tests  for  linear  behavior  for  various  q.  Second,  the  partition  function 
r(q)  is  computed  as  the  slopes  of  a  least  square  linear  fit  of  log  S^n\  Finally,  the  Legendre 
transform  r*(a)  of  r(q)  is  computed  following  (4.7).  Indicated  with  dashes  in  the  plots  of 
r(q)  and  r*(a)  of  pi,  are  the  corresponding  function  for  Mb,  providing  empirical  evidence  for 
(6.3),  (6.4),  and  (6.10). 

Corollary  6.3. 

<1]  =  E$,M  UAa)  =  U',M{a  +  1),  T^(a)  =  T^,M(a  +  1).  (6.11) 

This  is  impressively  demonstrated  in  Figure  7.  We  should  note  that  ip'  has  one  more 
vanishing  moment  than  ip  which  is  easily  seen  by  integrating  by  parts.  Thus,  it  is  natural 
to  analyze  the  integral  of  a  process,  here  the  distribution  function  M  of  the  measure 
p,  using  ip1  since  the  degree  of  the  Taylor  polynomials  of  M  are  typically  by  1  larger 
than  the  ones  of  p,  and  the  analyzing  wavelet  should  be  blind  to  these  polynomials. 

In  the  light  of  the  fact  that  (6.11)  holds  even  in  the  presence  of  chirps,  the  lesson 
to  retain  is  that  the  multifractal  analysis  of  p  =  M.'  through  ip  is  equivalent  to  the 
analysis  of  M  through  ip1  and  can  not  serve  to  discover  chirps.  It  is  notable,  then,  that 
corollary  6.2  uses  the  same  wavelet  for  both,  pb  and  Mb- 

Remark  6.4.  (Visibility  of  singularities  and  regularity  of  the  wavelet) 
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It  is  remarkable  that  the  Haar  wavelet  provides  us  with  the  full  spectra  of  both, 
the  binomial  M.h  and  its  distributional  derivative  pb-  This  seems  to  contradict  rules 
of  thumb  saying  that  a  wavelet  cannot  detect  degrees  of  regularity  larger  than  the 
number  of  vanishing  moments.  In  [78]  it  is  even  advocated  that  the  partition  function 
r(q)  should  be  used  only  for  q  G  [0, 1]. 

To  resolve  the  apparent  paradox  recall  the  peculiar  property  of  multiplicative  mea¬ 
sures  which  is  to  have  constant  Taylor  polynomials  in  the  points  of  interest.  So,  cascades 
will  reveal  their  scaling  structure  to  any  analyzing  wavelet  with  J  ip  —  0.  No  higher 
regularity,  i.e.,  f  tkip(t)dt  =  0  is  required  despite  the  fact  that  points  of  regularity 
a(t)  >  1  are  present.  For  another  example  of  a  signal  which  is  more  regular  than  the 
basis  elements  it  is  composed  of  see  [55]. 

The  mentioned  rule  has,  then,  to  be  corrected  to  read  as:  wavelet  cannot  detect 
singularities  in  points  with  approximating  polynomial  Pt  of  degree  higher  than  its 
regularity.  & 


7  Self-similar  Processes 

So  far  we  have  studied  multiplicatively  constructed  processes  which  exhibit  particular 
rescaling  properties  such  as  the  tree  structure  of  the  binomial  cascade.  Let  us  turn  now 
to  the  related  concept  of  so-called  ‘statistical  self-similarity’. 

7.1  Statistical  Self-similarity 

For  details  on  the  following  short  introduction  to  self-similar  processes  see  [100]  or  [97, 
p  309].  Here  is  some  notation  used  throughout:  /  is  either  IR,  1R+  =  {t  >  0}  or  {t  >  0}; 
V  is  a  self-similar  process  with  increments  U. 

Definition  7.1.  A  real-valued  stochastic  process  {Vt  :  t  G  1}  is  called  self-similar  with 
index  H  >  0  if  and  only  if  for  all  a  >  0 

V(at)  =  aHV(t),  (7.1) 

in  the  sense  of  finite  dimensional  distributions,  in  other  words,  if  and  only  if  for  any 
a  >  0,  any  m  G  IN  and  any  t\  <  t2  <,...,<  tm 

(V(ati), ,...,,  V(atm))  =  (aHV(h), ,...,,  aHV(tm )). 

{R(f)}  is  said  to  have  stationary  increments  if  and  only  if 

V(s  +  t )  -  V(t)  =  V(s)  -  y(o). 

A  complex  valued  process  is  self-similar  if  and  only  if  both,  real  and  imaginary  part  are. 
If  a  process  is  both  self-similar  and  has  stationary  increments  we  say  it  is  H- sssi. 
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It  follows  from  the  definition  that  V(0)  =  0  a.s.  for  any  self-similar  process.  Self¬ 
similar  processes  themselves  are  not  stationary;  as  a  matter  of  fact,  V(t)  is  self-similar 
if  and  only  if  e~tHV(et)  is  stationary.  It  is  helpful,  therefore,  to  add  the  assumption 
of  stationary  increments.  As  an  immediate  consequence  of  the  definitions,  an  H- sssi 
process  is  centered,  and  so  are  its  increments,  unless  it  has  infinite  expectation: 

IE R(l)  =  IE[V (1)  -  R(0)]  =  Wj\V (2)  -  y(l)]  =  (2h  -  1)IEI/(1). 

7.2  Examples 

If  one  assumes  finite  variance,  then  the  covariance  is  easily  computed  using  (7.1)  and 
stationarity  of  increments  as  follows 

r(s,  t )  =  1EV(s)V(t)  =  ^  (IE V{tf  +  IE V{sf  -  IE(H(t)  -  H(s))2) 

=  ^(\t\2H +  \s\2H -\t- s\2H)JEV(1)2.  (7.2) 

Since  this  expression  must  be  non- negative  definite  it  follows  that  0  <  H  <  1.  For 
H  =  1  one  Ends 

TE{V(t)  -  tV (l))2  =  WjV {t)2  -  2 1  ■  IEE(t)IEI/(l)  +  t2lEI/(l)2  =  0 

and,  thus,  V{t)  =  tV{l)  a.s.  for  all  t.  Often  one  studies  the  increments  of  an  H- sssi 
process.  Therefore,  let  us  set 


Uk  :=V(k  +  l)-V{k). 

Here,  we  choose  the  fixed  lag  or  step  1  for  convenience  only.  We  End,  using  the  expres¬ 
sion  for  r(s,  t) 

ru(k)  :=  IEo [U(n  +  k)U(n)]  =  -  (| k  -  1|2"  -  2\k\211  +  | k  +  1\2H )  .  (7.3) 

Gaussian  case 

As  mentioned  before,  the  auto-correlation  function  (7.2)  has  to  be  non-negative  definite 
which  restricts  H  to  (0, 1].  On  the  other  hand,  for  every  such  function  r  there  is  a  zero 
mean  Gaussian  process  defined  through  its  finite  dimensional  distributions  [97,  p  318]: 

IE  [GXp  X  ^  ^  6Xp  f  ~  EE  r(tk,ti)9k0i 

k= 1  \  k=  1  1=1 

It  can  be  shown  that  this  process  is  unique  and  that  it  is  always  possible  to  choose 
a  version  with  almost  surely  continuous  paths. 

Definition  7.2.  The  unique  Gaussian  and  H- sssi  process  with  0  <  H  <  1  is  called 
fractional  Brownian  motion  (fBm).  We  will  denote  it  by  Bh-  Its  increment  processes 
{B(t  +  5)  —  B(t)  :  t  e  1}  are  called  fractional  Gaussian  noise  (fGn). 
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The  case  H  =  1/2  is  notable  in  this  context  as  it  is  the  only  case  where  the  increments 
are  independent.  {Bi/2(t)}  is  called  Brownian  motion,  or,  to  make  a  clear  distinction, 
Wiener  motion  (WM).  FBm  has  been  introduced  by  Kolmogorov  and  studied  by  Levy, 
and  Mandelbrot  and  van  Ness  [71]  who  gave  the  following  representation  as  a  stochastic 
integral  over  Wiener  motion 

B«{t)  =  1  T  r [(i  -  sf+-'/2  -  (-s)r1/2]dBi/2 (s),  (7.4) 

1  [B  +  j)  J- 00 

where  (x)+  =  max(x,0). 

Of  particular  interest  here  are  the  fractal  properties  of  the  paths  of  fBm.  It  is  well 
known  that  its  paths  assume  almost  surely  everywhere  the  exact  Holder  regularity  of 
degree  H: 

Theorem  7.3  (Adler).  Let  [a,b]  be  a  compact  interval. 

a)  [5,  Thm.  8.3.1.]  For  every  rj  >  0  there  is  an  almost  surely  finite,  positive  random 
variable  A  such  that  with  probability  one  for  every  s,t  e  [a,  b]  and  |s  — 1\  small  enough 

\BH(t)  -  Bh(s)\  <  Als-tf-" 

b)  [5,  (8.8.26)]  For  each  rj  >  0  and  small  enough  h  >  0 

sup{|-Bf/(f)  —  Bh{s) I  :  \s  —  t\  <  h}  >  KhH+n 
for  every  t  G  [a,  b]  and  every  K  <  00,  with  probability  one. 

Stable  case 

If  we  drop  the  assumption  of  finite  variance  we  fall  on  the  7-stable,  H- sssi  processes  as 
the  next  larger  class.  A  symmetrical  7-stable  random  variable  V,  in  short  V  ~  S^S^cr), 
is  defined  through  its  characteristic  function: 

IE[exp(i#K)]  =  exp(— cr7|d|7  +  i/iO)  for  0  <  7  <  2. 

It  follows,  that  the  sum  V  +  V'  of  independent  SjS  variables  is  again  SjS  with 
<?l+v  =  av  +  aV" 

In  general,  no  closed  form  for  the  distribution  function  is  known.  Special  known 
cases  include  the  Normal  distribution  (7  =  2)  and  the  Cauchy  distribution  (7  =  1).  It 
is  known,  though,  that  the  tails  decay  with  a  power-law  for  0  <  7  <  2,  i.e.  there  is  a 
constant  C7  depending  only  on  7  such  that  [97,  p.  16] 

lim  A 7P[K  >  A]  =  lim  A 7P[K  <  -A]  =  C1  ■  a7.  (7.5) 

A— >00  A— >00 

Consequently,  the  absolute  moments  of  V  exist  only  up  to  order  7.  For  q  >  0:  IE|  W | 9  < 
00  if  and  only  if  q  <  7.  For  7  =  2,  the  distribution  is  a  Gaussian  which  has  finite 
moments  of  all  positive  orders. 


58 


R.  H.  Riedi,  Multifractal  Processes 


Let  us  denote  H- sssi,  S'yS'-stable  processes  by  Lnrr  Applying  the  definitions  we  find 
LHn(t)  -  Lh„{s)  ~  S~/S(\t  -  s\11  ■  uv (i) ) 

Since  we  may  write  L(t)  —  L(s)  +  L(s)  —  L{u)  =  L(t )  —  L(u),  the  increments  will  be 
independent  if  and  only  if  H  =  l/y,  which  generalizes  the  Gaussian  case  7  =  2.  A 
H- sssi,  S'yS- stable  process  with  H  —  I/7  is  called  Levy  stable  motion  and  we  denote 
it  by  Lh  :=  LHyi/H. 

Finally,  it  can  be  shown  that  necessarily  H  <  max(l/7, 1).  Except  for  the  Gaussian 
case  7  =  2  the  H- sssi,  7-stable  motion  are  not  unique  and  their  paths  are  discontinuous 
with  probability  one.  For  H  <  1,  a  class  of  fractional  stable  motions,  so-called  linear 
fractional  stable  motions,  can  be  constructed  explicitly  through  integral  representations 
similar  to  (7.4),  where  the  exponents  of  the  integration  kernel  are  now  H  —  I/7  instead 
of  H  —  1/2  and  where  the  integrating  process  is  Levy  stable  motion  instead  of 
Wiener  Brownian  motion  B\i 2. 


7.3  Multifractal  Properties  of  Self-similar  Processes 

Deterministic  Envelopes 

For  an  H- sssi  process  V ( t )  we  find,  due  to  V (0)  =  0, 

2-n»in)  _  | V({k  +  l)2-n)  -  V (k2~n)\  =  \V(2~n)\  =  2~nH\V(l)\  =  2-nH2~a°0> .  (7.6) 


More  generally: 

Lemma  7.4.  For  an  H-sssi  process,  the  singularity  exponents  07” \  h^\  and 
the  scaling  law 

s(n)  A  pj  +  (l/n)so0)  (s  G  {a,  h,  u>}). 

Furthermore,  IEn[2-9So is  finite  exactly  for  q  <  q  <q,  where  for  a(fl>  and 


obey 

(7.7) 


(?,  o) 


and  for  h 


(n) 

k 


(-1,00)  for  fBm, 

(—1,7)  for  S'yS  motion, 


q  =  cx)  for  fBm, 

q  <  7  for  S'fS  motion, 

q  =  —00  for  Levy  stable  and  Wiener  motion. 


(7.8) 


(7.9) 


Proof 

Property  (7.7)  follows  for  afi1  from  (7.6),  and  similarly  for  lifi}  using  continuity  of 
paths.  For  v/fl>  see  [24],  Property  (7.8)  is  well-known  (compare  (7.5))  since  2~nak  and 

(n) 

2 ~nwk  are  the  moduli  of  Gaussian,  resp.  stable  variables  [24], 

Next,  let  us  discuss  q  for  lif  '1 .  For  SyS-stable  motion  the  obvious  estimate  2~nak  < 

,  (n) 

2~nnk  implies  immediately  that  q  <  7.  For  fBm  we  note  that  [59,  1cm.  12.2.1]  actually 
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shows  that  the  tail  probability  P[supfer01i  \Bn(t)\  >  x]  is  smaller  than  4  exp  (—ex2), 
which  implies  q  =  oo. 

Finally,  still  for  h^\  q  =  —  oo  can  be  derived  at  least  for  Levy  stable  motion  Lh 
and  Wiener  motion  B\j 2  =  L2  using  stationarity  and  independence  of  increments  as 
follows: 


P[sup  \LH(t)\  <  e]  <  P[\Lh(*±L)  -  Lh(£)\  <  2s  for  k  =  0, .  .r3 m  -  1] 
o<«<i 

=  P[\LH(l/m)\  <  2e]m  <  (cm11  ■  e)m 

for  some  constant  c  and  all  m  (for  the  Gaussian  case  see  also  [48]).  This  is  sufficient 
to  show  that  all  moments  of  negative  order  exist.  We  believe  that  q  =  —  oo  also  for 
fBm.  0 

Using  (7.7)  one  finds  immediately 


2n  — 1 


En  [S(nW)j  =En^2 


(n) 

-nqsl  _  2r 


2-nHqJE< 


k= 0 


(0) 
-Is  0 


and,  with  (7.8)  and  (7.9) 


T,(q)  := 


lim 

n^oo 


logg Eng%) 

—n 


qH  —  1  for  q  <  q  <  q 
— oo  else. 


(7.10) 


For  the  computation  of  F(a)  let  us  denote  the  density  function^  of  the  distribution 


„(0) 


of  2  so  by  <hs.  The  P-sssi  property  (see  lemma  7.4)  and  the  mean  value  theorem  yield 


Pn[|4n)  ~a\<e]  =  Pn[n(a  -  H  -  e)  <  s^’ <  n(a  -  H  +  e)} 

2 -n(a-H-e) 

=  /  <hs(a;)da; 

J2~ri(a-H+s) 

__  2~n(a~H~£) —  2_2ne)<|>  (2~n^Xa’n~H^ 


M 


(7.11) 


where  xa,n  G  [a  —  £,  a  +  e] .  It  is  instructive  to  compare  this  to  the  rescaling  property  of 
the  marginals  of  the  log-normal  binomial  cascade  (5.32). 

Applying  lemma  7.4  ^(rr)  behaves  as  x ~^1  at  zero,  resp.  as  x~^~l  at  oo.  From 
this  the  asymptotic  behavior  of  (7.11)  follows  immediately  provided  we  know  whether 
Xatn  —  H  is  eventually  positive  or  negative.  This  observation  resolves  the  cases  a  —  e>H 
and  a  +  e  <  H .  To  cover  the  case  a  =  H  let  us  note  that  trivially  P[2-n£  <  |V(1)|  < 
2n£]  — >  1.  In  summary, 

..  (  (a  —  H)q  +  o(e,  n)  if  a  >  H  +  e 

-  log2  Pn[|4n)  -  a|  <  e]  —  <  o(e,  n)  if  a  =  H  (7. 12) 

n  [  (a  —  H)q  +  o(e,  n)  if  a  <  H  —  e. 


^Recall  that  2  ao  *  =  |U(1)|,  2  ho  =  sup[_lj2]  | V (t) |  etc. 
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where  o(e,n)  <  (1  jn)  log2(l  — 2 ~2n£)  +  c-e  for  sufficiently  large  n.  In  the  case  a  <  H  —  e, 
e.g.,  we  may  choose  c  =  2q  +  4.  From  this,  a  formula  for  F(a)  follows  immediately. 
Moreover,  (7.12)  implies  also  F(a )  =  F(a)  since  limsup^^  and  liminf^oo  may  differ 
at  the  most  by  ce.  fnhnite  values  for  the  partition  function  and  the  straight  parts  of 
the  resulting  Legendre  transform  where  mentioned  in  [66,  77]. 

Corollary  7.5.  For  an  H-sssi  process  the  deterministic  partition  function  obtained 
from  the  coarse  singularity  exponents  a^\  h ^  or  is  given  by  (see  Figure  8) 

f  1  +  q(a  —  H )  for  a  <  H 

T*{a)  =  l  1  for  a  =  H,  (7.13) 

^  1  +  q(a  —  H )  for  a  >  H. 

Thereby,  q  and  q  are  provided  in  lemma  l.f.  Moreover,  F(a )  =  £(a)  (compare  (3.24)) 
and 

Fa  =  Fw  =  T*  =  T* ,  Fh  =  T*.  (7.14) 


dim(^) 

T*(a)  =  T*(a) 

1  - 

s\ope=VH/*  1  - 

1  - 

1 

-  slope=  l/Hy\slope=- 1 

Levy-stable 

H 

1 

H 

1  ' 

H 

1  - 

1  - 

\  1  - 

V 

' \slope=-1 

fBm 

H 

H 

1  ' 

H 

FIGURE  8.  The  theoretical  spectra  of  some  self-similar  processes.  Top  row:  1/H-stable  Levy 
motion  (H-sssi)  bottom  row:  fractional  Brownian  motion  Bu  (H-sssi) 

Dimension  based  Spectra 

For  fBm  a  result  by  Adler  [5]  (compare  theorem  7.3)  states  that  with  probability  one 
h(t)  =  H  for  all  t.  So,  for  a  <  H  we  have  almost  sure  equality  of  the  spectra: 

fBm:  dim(4“')  =  fh(a)  =  rft*(o)  =  Tk*(a)  =  |  ^  “  =  (7.15) 

In  this  light,  it  is  appropriate  to  call  fBm  mono-fractal.  For  the  interested  reader  we 
mention  that  it  is  possible  to  find  multifractal  structure  in  the  occupation  measure  of 
Brownian  and  stable  motion,  however  not  in  terms  of  pointwise  exponents  but  pointwise 
densities  (see  [25]  and  references  therein). 
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For  Levy  stable  motion  Lh  {H- sssi,  7  =  1/ //-stable  process,  7  <  2)  Jaffard  citeJ5 
computed  the  dimension  spectrum  (compare  Figure  8  and  (7.17)): 

Levy  motion:  dim (Z/j^)  =  |  °  <  (7.16) 

Grain  based  Spectra 

Wiener  Brownian  motion  Z/1/2  as  well  as  Levy  stable  motion  (//  =  I/7)  have  indepen¬ 
dent  increments.  Thus,  corollary  7.5,  theorem  3.14  and  corollary  4.1  combine  to  yield 
almost  surely 

Levy  motion:  //(a)  =  fs(a)  =  Fs(a )  =  r*(a)  =  T*(a )  (7.17) 

for  all  a  for  which  T*(o)  >  0  or  =  —00.  This  is  valid  for  the  singularity  exponents  cy'fi'1 , 
h!f  as  well  as  for  u/fi1  for  compactly  supported  wavelets;  (7.13)  provides  a  formula  of 
T*.  Counting  on  the  strong  decorrelation  of  wavelet  coefficients  of  ffim  [35,  57]  and 
keeping  in  mind  remark  3.15  one  may  hope  to  find  the  same  result  true  for  fBm  with 
H  1/2.  Numerical  evidence  supporting  this  point  is  found  in  [41], 

Partition  Function 

Due  to  the  linearity  of  T(q)  we  are  able  to  compute  r(g)  for  any  H- sssi  process  assuming 
only  that  T(q )  is  finite  for  some  negative  and  some  positive  q.  We  note  first  that  for 
each  n,  S'^n^(0,ca)  counts  the  number  of  non- vanishing  increments  which  is  2n  almost 
surely.  Thus,  r(0)  =  —1  almost  surely.  Similarly,  T( 0)  =  —1  always.  Now  observe  that 
r(g)  must  be  concave,  yet  lie  ‘above’  T(q)  due  to  lemma  3.9.  This  is  only  possible  if 
r(g)  is  linear  as  well.  In  conclusion,  if  q  <  0  <  q,  then  with  probability  one: 

//-sssi:  T(q,oj)  =  qF[  —  1  for  all  q  <  q  <q.  (7-18) 

Almost  sure  bounds  can  be  derived  for  other  q  values  using  that  r(g)  is  concave  and 
that  r*(a)  is  positive  (compare  [62]). 

Failure  of  the  multifractal  formalism? 

A  stark  discrepancy  in  the  decreasing  part  between  the  local  spectrum  based  on  Holder 
exponents  dim(ZZj/)  (often  considered  the  ‘true  spectrum’)  and  the  deterministic  en¬ 
velope  T*w  based  on  increment  or  wavelet  exponents  afi'1  or  u/fi'1 ,  respectively,  is  im¬ 
mediately  apparent  from  Figure  8.  One  could  call  this  a  failure  of  Tfw  to  spot  the  true 

Hdlder  continuity.  Technically  speaking,  the  reason  for  this  feature  is  the  fact  ^  and 
v/fi'1  are  centered  (zero-mean)  random  variables  of  cn,  actually  Gaussian,  resp.  stable 
variables  with  considerable  probability  around  zero. 

Keeping  in  mind  (7.17)  and  (7.18),  however,  one  has  to  acknowledge  that  Tfw  nev¬ 
ertheless  contain  crucial  scaling  information  which  is  relevant  in  real  world  application. 
The  part  of  these  spectra  forming  the  line  of  slope  —1  represent  the  considerably  high 
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chance  to  observe  large  a ^  and  w ^  at  resolution  2_n,  i.e.,  increments  and  wavelet 
coefficients  of  magnitude  smaller  than  2~nH . 

In  order  to  obtain  information  on  the  true  Holder  continuity  from  wavelet  exponents 
w one  has  to  employ  the  lines  of  maxima  (see  [78,  79,  7,  23]).  Taking  local  maxima  of 
the  modulus  of  wavelet  coefficients  one  avoids  exceptionally  small  ;  the  correspond¬ 
ing  partition  function  comes  then  closer  to  Th  since  less  mass  in  concentrated  around 
zero,  allowing  for  more  moments  of  negative  order  to  be  finite.  The  interested  reader 
might  want  to  compare  this  remedy  with  the  use  of  which  avoids  exceptionally 
small  increments  for  measures  on  fractal  support  (5.36). 

In  conclusion,  one  should  keep  in  mind  this  artifact  of  multifractal  analysis  based  on 
data  such  as  increments  which  is  of  zero  mean.  Though  some  information  on  the  local 
singularity  structure  of  the  process  might  be  hidden  to  such  an  analysis,  the  resulting 
T*  nevertheless  corresponds  to  true  scaling  behavior  of  the  process  in  the  sense  of  (7.17) 
which  might  be  relevant.  In  this  context,  it  becomes  essential  to  detect  divergence  of 
moments  from  finite  data,  an  issue  which  is  under  investigation  at  present  [39]. 


7.4  Long-Range  Dependence  (LRD) 

The  rigid  correlation  structure  of  fGn  (7.3)  is  somewhat  restrictive  for  modelling  pur¬ 
poses.  However,  the  fact  that  its  auto-correlation  decays  very  slowly  has  been  found 
to  be  an  important  feature  in  itself  and  has  inspired  weaker  notions  of  “similarity  on 
all  scales”  in  terms  of  second-order  statistics  only.  The  novice  reader  should  avoid  a 
common  confusion  with  these  weaker  notions  and  be  aware  that  they  are  defined  in 
terms  of  the  increment  processes  on  varying  lags  rather  than  in  terms  of  the  process  Y 
itself. 


Second  Order  Scaling 

It  is  easy  to  see  that  (7.3)  decays  like  rjj(k)  ~  k2H~2.  For  1/2  <  H  <  1,  the  correlation  is 
strictly  positive  and  decays  so  slowly  that  it  is  not  sumable  over  k.  A  process  X  -  which 
we  think  of  as  being  the  increment  process  of  Y  -  with  this  property  (ff2k  TX{k)  =  oo) 
is  said  to  exhibit  long  range  dependence  (LRD),  since  it  exhibits  strong  correlations 
at  large  lags.  LRD  can  be  equivalently  characterized  in  terms  of  the  behavior  of  the 
aggregated  processes 

(fc+l)m— i 

X{m\k)  :=  -  V  X(i).  (7.19) 

m 

i=km 

As  [19]  shows,  rx[k )  ~  k2H ~2  is  equivalent  to  var(A"(/c))  ~  rn2^2HVcir(X(m'> [n] ) .  So, 
we  are  lead  to  define: 

Definition  7.6.  We  say  that  X  is  asymptotically  second  order  self-similar  with  pa¬ 
rameter  Hvar  if  the  following  limit  exists: 


Hv ar  1  + 


1  ,  logfvaHX^VvaHAG1))] 

-  lim  - - - 

2  m— >oo  log  m 


(7.20) 
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As  an  example  let  us  mention  the  increments  of  an  H- sssi  process  V  (of  finite  vari¬ 
ance):  setting  X(k)  =  U{k)  =  V(k  +  1)  —  V(k)  one  finds  Hva r,x  =  H  due  to  (7.1)  or 
(7.6)  (compare  (1.17)). 

Alternatively,  LRD  can  be  measured  through  spectral  properties,  since  rxik )  ~ 
k2H~2  is  equivalent  to  a  power-law  behavior  of  its  Fourier  transform,  the  power  spec¬ 
trum,  i.e.,  Yx{y)  =  tF(rx)(v)  —  u~ly2H^1'1 .  The  spectral  estimation  using  the  wavelet 
transform  of  X  developed  in  [4]  translates  this  into 

var {Cfk)  oc  2~j{2H-1l  (7.21) 


One  may  adopt  the  following  terminology: 

Definition  7.7.  X  is  asymptotically  wavelet  self-similar  with  parameter  Hwave  if  the 
following  limit  exists: 


H , 


-  lirn  -  log2  var((7j  fc). 
Z  j->o o  J 


(7.22) 


Also,  Hwave  can  be  considered  an  estimator  of  H  for  H- sssi  processes.  Similar  as 
before,  Hwave  =  H  for  fGn  [35,  57].  It  should  be  noted  here  that  the  wavelet  coeffi¬ 
cients  of  fGn  are  typically  much  less  correlated  than  those  of  the  underlying  sampled 
fBm  process.  Thus,  wavelet-based  estimators  for  H  compare  favorably  with  standard 
estimation  techniques  [4]  and  is  particularly  superior  to  FTvar  [99]. 


Multifractal  Envelope  and  LRD 

The  multifractal  scaling  exponent  T( 2)  of  a  process  Y  is  closely  related  to  the  LRD 
parameter  H ,  since  both  measure  the  power-law  behavior  of  second-order  statistics. 
More  precisely,  T( 2)  captures  the  scaling  behavior  of  the  second  sample  moments, 
while  H  captures  the  decay  of  the  covariances. 

From  Section  1.5  we  recall: 


Lemma  7.8.  For  a  process  Y  with  zero-mean  increments 


1 1  va  r 


Ta,y(  2)  +  l 


(7.23) 


For  fBm,  this  is  in  agreement  with  (7.13). 

Multifractal  processes  defined  through  multiplication  on  a  tree  such  as  the  bino¬ 
mial  X{k)  =  are  not  second-order  stationary.  Hence,  LRD  cannot  be  defined 

through  the  decay  of  the  auto-covariances.  However,  alternative  fractal  properties,  such 
as  the  decay  of  aggregate  variances  (7.20)  or  wavelet  coefficients  (7.22)  —  which  are 
equivalent  to  LRD  in  the  presence  of  second-order  stationarity  —  can  still  be  defined 
and  calculated. 

As  a  further  difficulty,  processes  obtained  from  cascades  have  positive  increments 
X(k),  so  that  the  argument  (1.18)  of  Section  1.5  using  the  variances  has  to  be  corrected 
to  read 


2-2i2(i-n)(i+T(2))  ^  E[|xM|2]  =  var(x(m))  +  IE[X(m)]2  ~  var(X)2i(2^~2)  +  1E[X]2, 
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noting  that  IE^Y^A]  is  independent  of  the  scale  m.  Since  2 H  —  2  <  0  we  may,  thus, 
still  expect  the  same  relation  (7.23),  at  least  in  the  limit  of  very  fine  resolution  (small 
m  and  i). 

The  variance-time  plot  method  (7.20)  is  known  to  be  an  unreliable  (but  simple)  esti¬ 
mator  of  LRD  behavior  [99],  while  the  wavelet  method  (7.22)  is  more  robust.  Recalling 
that  TWtx{ 2)  is  defined  in  terms  of  wavelet  energy,  i.e.,  the  variance  of  the  wavelet  co¬ 
efficients  of  X,  we  find  by  stationarity  var (Cjk)  =  IE[2J'|(7j  fc|2]  ~  2_h1+T(2))  — 

2-j(2+T(2))j 

=  rmJf(22)  +  3.  (7.24) 

Example  7.1.  (LRD  parameter  of  Binomial  Cascades)  Together  with  (6.3)  we 
hnd,  e.g.,  for  the  binomial  Y  =  Xib: 


H, 


wa,ve,A4b 


Tw,nb(2)  +  3  _  TMb( 2)  +  1 


=  if. 


var,A^t  • 


(7.25) 


* 

This  formula  has  interesting  similarities  with  (6.2). 

Example  7.2.  (When  scaling  is  only  asymptotic)  The  examples  presented  so  far 
in  this  section  are  processes  with  ‘perfect’  or  ‘omni-scaling’  such  as  fBm  and  cascades, 
meaning  that  the  scaling  parameters  H  and  T(q)  can  be  obtained  from  comparing 
statistics  of  any  two  scales.  In  real  world  applications  one  often  requires  more  flexibility 
and  matches  the  LRD,  resp.  multifractal  parameters  only  in  the  limit  of  large,  resp. 
small  scales.  Potentially,  the  two  matching  procedures  can  be  performed  independent 
of  each  other. 

Linear  processes  such  as  FARIMA  (or  ARFIMA)  [100]  allow,  e.g.,  to  match  short 
term  correlations  in  addition  to  LRD.  As  a  first  step  towards  decoupling  LRD  and 
multifractal  scaling  we  mention  here  cascades  with  scale  dependent  multipliers.  Such 
processes  were  used  in  [89]  to  model  internet  traffic.  In  doing  so  we  also  demonstrate 
how  misleading  time-domain  estimators  such  as  (T(2)  +  l)/2  can  be  when  only  few 
scales  of  resolutions  are  available  and  a  large  mean  ‘masks’  the  fine  scale  variability. 
Let  us  add  that  ‘fBm  in  multifractal  time’  which  we  introduce  in  the  remaining  section 
provides  an  attractive  alternative  to  model  LRD  separately  from  fine  scale  properties. 

Consider  a  binomial  cascade  supported  on  [0, 1]  with  (i)  and  (ii),  but  where  instead 
of  (iii)  the  distributions  of  M jcn'>  depend  on  the  scale  n  as  follows: 


var(MjJ")) 


22-2gvar(M^~1}) 
4var(il4n_1))  +  1  ' 


(7.26) 


This  is  arranged  such  that  the  Haar  wavelet  coefficients  scale  exactly  as  in  (7.21).  This 
powerlaw  can  be  easily  derived  using  (5.13)  (see  below  for  a  rough  argument  or  see 
[89]).  In  short,  it  is  fair  to  say  that  this  process  has  a  scaling  of  the  LRD  type  with 
exponent  H . 
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To  fully  determine  the  distributions  assume  that  M^n  1}  follows  a  symmetrical  (3- 
distribution  on  [0, 1]  (compare  Example  5.1  or  [89]).  To  start  the  iteration,  we  assume 
a  value  of  E[(M^)2]  ~  1/4  (but  different  from  1/4)  or  var(M^C)  ~  0.  This  reflects 
typical  Endings  with  high  resolution  network  traces  where  there  is  often  little  variation 
in  the  data  on  high  aggregation  levels  in  the  order  of  several  minutes. 

Iterating  (7.26)  we  find  for  the  first  few  n  that  IE [(M^C)2]  ~  1/4  and  var (M^)  ~ 
a2(22~2H)n,  which  is  small  for  small  n.  Only  for  very  fine  scales  (large  n)  will  1E[(M^)2]  = 
var (M^1  )  +  1/4  converge  to  the  limiting  2~2H .  All  this  will  affect  the  estimation  of  the 
LRD  parameter  through  (T( 2)  +  l)/2  from  limited  data  as  we  explain  now. 

Motivated  by  (5.24)  we  use  —  (1/n)  log2(2n  YYi=i  E[(M^)2])  as  an  estimator  for  T( 2). 

If  measurements  are  available  only  for  small  n  this  will  result  in  T( 2)  ~  —(1/n)  log2(2n4~n 
1  since  convergence  is  slow,  and  (T( 2)  +  l)/2  ~  1.  Only  in  the  limit  n  oo  will 
one  see  T( 2)  ~  —  (1/n)  log2(2"2~2J?n)  =  2 H  —  1  and  discover  the  ‘correct’  value 
(T( 2)  +  l)/2  ~  H. 

This  is  in  stark  contrast  to  the  ability  of  the  wavelet  based  estimator  Hwstve  to  pick  up 
H  even  from  observing  only  few  scales  (recall  that  (7.21)  holds  exactly).  Indeed,  we  ar¬ 
gue  that  (7.21)  should  hold  even  when  n  is  still  small  (compare  [89]  for  a  rigorous  proof 
of  the  exact  powerlaw).  According  to  (5.13)  and  using  M^1  ’  —  =  2M2]TK)  —  1 

we  may  translate  the  above  approximations  into 

IE [Cltk\  =  2"1E[(M^1))2]  •  ■  ■  E[(M^n))2]  •  4var(M^n+1))  ~  2n4“n4(22-2H)n  =  A(21~2H)n 

(7.27) 

as  claimed.  Consequently,  the  wavelet  energy  reveals  Hw ave  =  H  even  for  small  scales. 

We  see  the  reason  for  this  phenomenon  in  the  fact  that  T( 2)  is  based  on  non-centered 
2nd  moments:  for  this  estimator  of  LRD  the  scaling  behavior  is  at  coarser  scales  hidden 
behind  the  mean  of  the  process.  In  further  support  of  this  conclusion  we  mention  that 
also  ifvar  ~  H  even  when  estimated  from  few  resolutions  (n  small):  indeed,  the  ‘first 
order’  approximation  IE[(M^)2]  ~  1/4  {i  =  1, . . . ,  n  —  1),  E =  1/4+var (M^) 
results  in  var(Af(2'n))  =  E[(M^1})2]  •  •  •  E[(M^n))2]  -  (l/2n)2  ~  l/4n-1var (M^)  ~ 
u2/ 4  •  2~~2Hn,  which  is  the  claimed  scaling. 

As  a  final  remark  we  add  that  with  the  ‘initial  condition’  E[(M^)2]  =  2~m  scaling 
is  perfect  and  (7.25)  holds  again. 


8  Multifractal  Time  Change 

So  far,  we  have  given  examples  of  processes  which  are  either  statistically  self-similar 
and  ‘mono- fractal’  or  multiplicative  and  almost  surely  increasing^  and  multifractal.  It 
is  clear  that  there  is  need  to  marry  the  two  in  order  to  have  processes  with  richer 
structure.  To  quote  Mandelbrot  the  idea  towards  this  goal  is  ’’simple  but  inevitable 


I  Our  multipliers  are  positive.  Vector-valued  multipliers  can  be  treated  but  lead  to  technical 
difficulties  [31]. 
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after  the  fact” :  multifractal  time  change. 

Indeed,  ‘Brownian  motion  in  multifractal  time’  appeals  through  the  simplicity  in 
which  it  combines  attractive  features  such  as  Gaussianity,  LRD  and  multifractal  struc¬ 
ture.  Arguing  that  trading  occurs  in  ‘multifractal  time’,  [68,  69]  introduced  this  process 
as  a  model  for  stock  exchange.  Here,  we  are  exploring  all  its  multifractal  spectra,  as 
well  as  more  general  ‘multifractal  subordinations’.  Since  most  of  our  analysis  can  be 
done  pathwise,  let  us  start  with  a  deterministic  example.  Recall,  that  we  denote  the 
spectra  of  a  process  Y  computed  with  singularity  exponents  s ^  by  fS}y  etc.  when 
confusion  is  possible. 

8.1  Weierstrass  Function  in  Multifractal  Time 

The  Weierstrass  function 

W(t)  :=  X~kH  sm(Afct)  (8.1) 

has  originally  been  introduced  as  an  example  of  a  nowhere  differentiable  function.  In  its 
randomized  version  it  was  proposed  by  Mandelbrot  as  an  efficient  and  simple  model  of 
fBm.  Similarly  to  fBm,  W ( t )  is  a  mono- fractal,  more  precisely,  its  local  Holder  exponent 
hw(t )  is  everywhere  equal  to  H .  More  pointedly,  there  are  c  and  d  such  that 

\W(t  +  5)-W(t)\<c8H  \W(t  +  e)  -  W{t)\  >  deH  (8.2) 

for  all  t,  all  6  and  some  e  <  5. 

In  order  to  leave  uniform  Holder  continuity  behind  and  moving  towards  multifractal 
scaling  we  proceed  as  follows: 

Definition  8.1.  Let  M.(t)  =  /x([0,t])  be  a  multifractal  distribution  function.  Then,  the 
Weierstrass  function  in  multifractal  time  W(MF)  is  given  by 

W(t)  :=  W(M{t))  (8.3) 

We  claim  that  W(MF)  has  a  significant  multifractal  structure.  Indeed,  assuming  that 
M.  is  continuous  in  a  neighborhood  of  t  we  find  a  5o(t)  >  0  such  that  for  5  <  8o(t) 

sup  |W(s/)  —  W(s)|  = 

[t-S,t-\-S] 

sup{  \W{u')  -  W{u)\  :  inf  M(s)  <  u  <  v!  <  sup  Al(s)} 

[t-6,t+S] 

<  c  sup  \M (s')  —  M(s)\u  (8.4) 

[t— 5, £+5] 

and  the  reverse  inequality  with  d.  From  this,  it  follows  immediately  that  for  all  e  >  0 
and  n  >  no(e,  t) 

\Y}m-H-h^(M)\<e.  (8.5) 


This  yields 
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Lemma  8.2  (Warp  formula  for  Weierstrass).  Assume  that  AA  is  piecewise  con¬ 
tinuous  on  [0, 1]  with  at  most  a  finite  number  of  discontinuities.  Then,  in  points  of 
continuity  of  AA 

hw(t)  =  H  ■  hM(t)  (8.6) 

and 

dim(^Sv)  =  di  m(Eh,M])-  (8-7) 

Moreover,  if  Ad  is  continuous,  then  it  is  actually  uniformly  continuous  on  [0, 1]  and 
it  is  easy  to  see  that  <50  and  n0  of  (8.4)  are  independent  of  t.  Thus,  we  may  relate  the 
grain  based  spectra  of  W  and  Ad  using  (8.5),  i.e.,  /yy(a)  =  /x(a/ih).  Their  partition 
functions  are  connected  by  a  similar  formula  as  we  are  about  to  develop.  Referring 
back  to  (8.4)  we  find,  e.g.,  for  positive  q,  and  n  >  no 


2"-l 


2n-l 


£  <  cqH  2~qnH-h<?)W  =  cqHS(hn)M(qH). 


k= 0 


k= 0 


Similar  estimates  using  the  reverse  of  (8.4)  show  that  Sjf^fq)  and  Sjf^^q)  are  equal 
up  to  a  factor  which  is  bounded  independently  of  n.  So,  if  tm  (qH)  assumes  a  limit  as 
in  (3.12),  then  so  does  Tyv(q).  In  summary 

Lemma  8.3.  Assume  that  Ad  is  continuous  on  [0, 1].  Then, 


Th,w(q)  =  Th,M(qH)  fh,w(a)  =  fh,M(a/H ).  (8.8) 

These  lemmas  say  nothing  more  than  what  was  already  apparent  in  (8.2):  the  Weier¬ 
strass  function  stretches  distances  locally  in  a  very  uniform  way  in  the  form  of  a  pow- 
erlaw  with  exponent  H.  In  W(MF)  in  particular,  it  stretches  the  multifractal  time 
increments,  multiplying  local  Holder  exponents  by  H .  As  a  consequence  we  have: 

Theorem  8.4  (Multifractal  formalism  for  Weierstrass).  Let  W  be  defined  as  in 
(8.3)  and  assume  that  AA  is  continuous.  Consider  the  singularity  exponents  h\f\  Then, 
for  every  a  with 

di  miEMH])  =  TXi(a/H ) 

equality  holds  in  the  multifractal  formalism: 

dim(£,{J)  =  fw(a)  =  r^(a)  =  r*M{a/H).  (8.9) 


More  delicate  questions  of  multifractal  analysis  find  again  a  simple  answer  provided 
AA  is  continuous.  E.g.,  if  E^rr  has  positive  finite  y-dimensional  Hausdorff  measure 
then  Eyj  has  positive  finite  7- dimensional  Hausdorff  measure  where  7  =  dim = 
dim(ifj/y ),  etc.  Similarly,  it  follows  that  W  has  a  rich  multifractal  structure  if  Ad  has 
one.  By  this  we  mean  that  the  sets  Ey(\  are  dense  on  the  unit  interval.  A  smooth,  strictly 
concave  partition  function  Tyy  can  be  taken  as  an  indicator  of  highly  interwoven  sets 
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The  particular  form  of  (8.9)  and  (8.8)  allows  to  separate  the  Holder  continuity  H  of 
the  Weierstrass  function  from  the  multifractal  spectrum  of  the  ‘time  change’  M,  even 
numerically.  Indeed,  the  slope  of  the  tangent  of  the  W-spectra  through  the  origin  have 
slope  1  / H  since  ry^{l/H)  =  tm(1)  =  0.  Alternatively,  one  may  look  for  the  zero  of  ryy 
in  order  to  determine  H .  Once  this  is  accomplished  follows  easily. 


8.2  Self-similar  processes  in  multifractal  time 

Following  the  example  of  the  last  section  it  is  now  easy  to  define  processes  with  rich, 
ubiquitous  multifractal  structure,  yet  attractive  statistical  properties  such  as  zero  mean 
marginals  and  second  order  stationarity,  simply  by  multifractal  time  change.  In  order 
to  be  able  to  compute  multifractal  spectra,  one  needs  some  scaling  properties  of  the 
subordinated  process,  of  course.  To  this  end,  let  V  denote  a  H- sssi  process  and  M(t) 
an  almost  surely  increasing  process  from  which  we  assume  to  know  certain  of  its  mul¬ 
tifractal  spectra.  Then,  set 

V  :=  V(M(t)).  (8.10) 

As  always,  the  deterministic  partition  function  T(q)  is  the  most  simple  of  the  spectra 
to  calculate.  Indeed,  the  computation  for  H- sssi  processes  carried  out  in  the  introduc¬ 
tion  (see  (1.24))  can  easily  be  generalized. 

Theorem  8.5.  Let  V  denote  an  H-sssi  process  and  M  an  independent  time  change. 
Let  V  :=  V(M(-)).  Then,  with  q  and  q  as  defined  in  lemma  7.4, 


Tv(q) 


TM(qH)  if  q<  q  <  q 
— oo  else. 


(8.11) 


Here,  the  singularity  exponents  may  be  chosen  to  be  a ^  without  further  assumption. 
For  tifi'1  the  result  holds  provided  M  has  almost  surely  continuous  paths.  Finally,  in 
the  case  of  a  binomial  warp  time  Mb  as  in  Section  5,  (8.11)  holds  also  for  for  any 
admissible,  analyzing  wavelet  supported  on  [0, 1]  under  some  technical  condition. 

Note  that  Ty(g)  may  assume  the  value  —  oo  even  for  q  <  q  <  q  depending  on  whether 
Tj^(qH)  is  finite  or  not. 

Proof 

First,  assuming  continuous  warp  time  M  we  may  deal  with  suprema  as  in  (8.4).  This 
takes  care  of  the  case  h  f  } . 

Let  us  now  extend  (1.24)  to  v/f1'1 .  When  conditioning  on  M  we  obtain  by  H- sssi  and 
(5.16) 

m  -  v(k 2-")  \M  a  (m£>  ■  ■  ■  m<;>  ■  -  *.))"  ■ 
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The  admissibility  of  the  wavelet  (f  if  =  0)  allows  now  to  write 

2 n/2C„,tJV)  \M  =  2”  f  V'(2 "t  -  kn)V(t)dt  \M 

J  lb. 

Kn 

=  2n  J ^  if(2nt  -  kn)  (V(t)  -  V(kn2~n ))  dt  | ^ 

^kn 

=  (M<:>  ■  •  •  <)"**  (M?)  ■  V(l).  (8.12) 

Hereby,  An>kn(Ai^)  stands  for  the  integral 

1  2n  — 1 

f  if(t)(Mb{t))Hdt  =  lim  V  iP(k2-n)(Mb(k2~n))H  ■  2~n  (8.13) 

/n  71— >00  '  J 

fe=0 

and  is,  actually,  nothing  but  the  wavelet  coefficient  Co,o  of  {M.^1,kn^)H .  Since  A4^l'kn)  =Mb, 
we  have  En[|An)fcn  |9]  =  IEf)[|C'o)o(A/l^)|'3]  and  the  integral  term  will  not  affect  the  scal¬ 
ing  of  moments,  provided  IEn[|Coio(A/f^)|IJ]  is  finite.  This  is  one  part  of  our  technical 
assumption. 

Similarly,  JEn[|V(l)9]  will  not  affect  TvjV  as  long  as  it  is  finite.  This  is  taken  into 
account  by  restricting  q  to  (q,q). 

Comparing  with  the  derivation  (5.24)  one  finds  now  easily  that  Twy(q)  =  Ta)Mb(qH). 
Assuming  now  also  that  En[|Co,o(Al6)|9'f/]  is  finite,  which  is  the  second  part  of  the 
technical  assumption,  then  TaMb(qH )  =  Tw  Mb(qH)  by  (6.3).  So,  (8.11)  is  established. 

❖ 


8.3  FBm  in  multifractal  time 

The  most  obvious  example  of  a  ‘warped  motion’  is  fractional  Brownian  motion  in 
multifractal  time  FB(MF),  i.e.,  B  =  It  is  first  mentioned  in  [68,  69]  which 

pioneered  it’s  application  to  stock  markets.  FB(MF)  has  particularly  nice  properties. 
It  must  be  distinguished  from  the  ‘multifractional  Brownian  motion’  of  [85]  which  has 
a  ‘slowly  varying’  h(t),  non-ubiquitous  singularities,  a  piecewise  linear  T*,  and  strict 
inequality  in  the  multifractal  formalism  [62],  i.e.,  an  excellent  model  for  non-stationarity 
rather  than  for  multifractality. 

Pathwise  Spectra  of  warped  fBm 

For  any  process  V  which  has  strict  and  uniformly  Holder  continuous  paths  such  as 
fBm,  the  pathwise  spectra  are  as  easy  to  compute  as  for  the  multifractal  Weierstrass 
function.  Indeed,  applying  Adler’s  result  (see  Theorem  7.3)  for  r)  =  1/m  (■ m  G  IN) 
instead  of  (8.2)  we  find  that  with  probability  one,  for  all  m  G  IN  there  is  a  finite, 
positive  Am  and  n0(m )  G  IN  such  that 

“log  2<yAm)(H  -  1  /m)h^(M)  <  h^{B)  <  ^  log  2(Am)(H  +  1  (8.14) 
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for  all  t  and  all  n  >  ni(m,t).  As  usual,  =  kn(t )  =  [t2n\  as  in  definition  5.1.  Here, 
rii  has  to  be  taken  large  enough  such  that  M.  is  continuous  in  [{kn  —  l)2~n,  (kn  +  2)2"n] 
(compare  (8.4))  and  such  that  lA^s')  —  A4(s)|  <  2~n°  for  all  s,  s'  G  [{kn  —  l)2_n,  (kn  + 
2)2~n],  To  avoid  complications  we  should  assume  here  that  time  change  At  and  ffim 
Bh  are  independent. 

Assuming  that  A4  has  almost  surely  continuous  paths  this  allows  to  conclude  imme¬ 
diately  that  with  probability  one, 


hB{t)  —  H  • 


for  all  t  in  a  compact  interval.  When  M.  is  continuous,  however,  then  n \  can  be  chosen 
independently  of  t  with  probability  one,  and  we  may  compare  grain  based  spectra  and 
the  partition  functions  in  the  same  way  as  for  the  warped  Weierstrass  function  W(MF). 
So,  lemma  8.3  and  theorem  8.4  hold  for  warped  ffim  FB(MF)  mutatis  mutandis: 

Theorem  8.6  (Warp  formula  for  fBm).  Let  Bh  denote  fBm,  let  Jvl  be  an  almost 
surely  continuous  time  change  independent  of  Bh  and  set  B(t)  :=  Bni^Aft)).  Then, 
for  almost  every  path  (8.7)  and  (8.8)  hold  with  W  replaced  by  B.  In  particular,  for 
almost  every  path 


dim  (Fig1)  =  f  13(a)  =  Tg(a )  =  r*M(a/H)  (8.15) 


for  all  a  with 


dim  (E[^h])  =  T*M(a/H), 


where  all  these  spectra  are  taken  with  respect  to  the  Holder  exponents  h 


(n) 

k 


Estimation  through  wavelets 

While  theorem  8.6  is  an  appealing  theoretical  result,  it  is  quite  unaccessible  from  a 
numerical  point  of  view  because  suprema  of  continuous  time  processes  are  hard  to 
estimate  from  a  discrete  sampling.  It  comes  as  a  remedy  that  wavelets  can  be  used  for 
estimation  purposes  as  we  are  going  to  develop  in  a  special  case  now. 

Therefore  let  us  assume  that  the  warp  time  is  a  binomial  cascade  satisfying 
(i)-(iii)  (see  (5.5)).  Let  if  be  an  admissible,  analyzing  wavelet,  i.e.,  f  if  =  0,  supported 
on  [0, 1].  Exploiting  the  uniform  Holder  continuity  of  ffim  (see  (8.14)  or  Theorem  7.3) 
and  the  strong  scaling  property  of  M. &  (see  (5.14))  we  will  lay  open  the  multiplicative 
structure  of  the  wavelet  coefficients  Cn^n{B)  in  a  result  analogous  to  lemma  5.4.  This 
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is,  of  course,  much  stronger  than  (8.12)  which  dealt  only  with  distributions: 
2  n/2Cn,kn(B)  = 


Hereby,  we  used  ~  to  indicate  equality  up  to  a  bounded  error,  i.e.,  A  and  K  in  the 
notation  of  theorem  7.3.  Also,  Anjkn(A4fr )  is  as  in  (8.13)  and  we  assume  it  to  be 
uniformly  bounded  away  from  zero  for  H'  —  1  as  well  as  all  H'  close  to  H.  Taking  log2 
and  dividing  by  —n  we  find  the  analogue  to  (8.14),  i.e., 

«  H  ■  w^(Mh) 

up  to  an  error  which  will  become  uniformly  small  as  n  increases.  In  short,  wB{t)  = 
H  ■  w Mb (t),  ,/u',e(a)  =  fw,Mb(a/ H)  e^c-  Since  for  a  binomial  the  wavelet  based  spectra 
coincide  with  the  ones  based  on  a ^  and  almost  surely  we  can  summarize  as 
follows: 

Theorem  8.7  (Wavelet  analysis  of  warped  fBm).  Let  A4b  be  a  binomial  time 
change  satisfying  (i)-(iii)  (see  (5.5)),  independent  of  Bn-  Consider  an  admissible, 
compactly  supported  analyzing  wavelet.  Assume  that  Anjk(A4(f  )  are  uniformly  bounded 
away  from  zero  for  H'  =  1  and  H'  close  to  H .  Then,  for  almost  all  paths 

di  =  cMO  =  fwAa)  =  fhlB(a)  =  r*Aa)  =  ThAa )  =  TMb(a/H )  (8-16) 

for  all  a  with  (a)  >  0. 

This  is  impressively  demonstrated  in  Figure  3.  For  the  statistics  of  the  estimator  for 
TwAa 0  used  in  there  see  [40]. 

We  would  like  to  point  out  that  this  result  is  not  so  much  due  to  the  dyadic  structure 
of  both  the  binomial  warp  time  A4b  as  well  as  the  wavelet  coefficients  CU:k,  but  rather 
the  strong  rescaling  properties  of  M.b  and  fBm.  Using  sophisticated  tools  which  are 
now  standard  in  multifractal  analysis  it  is  very  well  possible  to  extend  this  result  to 
arbitrary  multi-nomial  warp  time  as  introduced  in  Section  5.6. 

LRD-parameter  for  warped  fBm 

An  immediate  consequence  of  (8.11)  is  a  formula  for  the  ‘LRD’  parameter  of  FB(MF). 
We  could  employ  (7.20)  or  (7.22).  However,  we  choose  to  calculate  correlations  exactly. 


r(") 
1k, i 


r(") 

1kn 


j(n) 

1kn 


ip (2 nt  -  kn)B(t)dt 

if(2nt  -  K)  (V{Mb{t))  -  V(Mb(kn2-n )))  d t 
if(2nt  -  K)  (. Mb(t )  -  Mb(kn2~n))H±v dt 


(<>  ■  ■  ■  K)  "M  [m  ( m to)  if 


(2 n/2Cnk  (. Mb))H±n  —  v — ~ — ttz —  * 
V  ,nX  ”  (A^u(MA\H±r) 


An,kn  (Mf^ 
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In  order  to  do  so,  we  develop  (B(t)  —  £>(s))2  and  take  expectations  as  in  (7.2).  First 
note  that 


IE[£>(t)2]  =  1EJE){BH{M{t))2\M{t)}  =  IE [M(t)2H]  ■  JF{BH(l)2}. 

Writing  a2  =  lEf-B#  (l)2]  for  short  we  obtain 

2 

IE [B(t)B(s)]  =  y  (jF[M(t)2H]  +IE[AI(s)2H]  -TE[\M(t)  -  Ad(s)|2H]).  (8.17) 

Let  us  consider  now  the  increment  process  for  a  fixed  lag  size  5  <  1,  i.e.,  fractional 
Gaussian  noise  in  multifractal  time  FG(MF): 

G(s)  :=  B(s  +  5)~  B(s)  =  Bh(M(s  +  5))  -  Bh(M(s)).  (8.18) 

Writing  =  M.{b)  —  M.{a )  for  short  we  find  from  (8.17) 

e  [aw  ■  aw  ]  =  yE  [iaWi2"  -  lO2"  -  iwf  "  +  IWWH  ■  (s.is) 

Lemma  8.8  (LRD  of  warped  fBm).  Assume  that  for  each  k  the  moment  of  order 
2 H  of  the  increments  A4(k8)  —  A4(mS)  depends  only  on  \k  —  m\  and  is  finite.  Then, 
the  discrete  time  FG(MF)  is  second  order  stationary  with  auto- correlation  function 


rg(k )  =  IEn[C7(n  +  k)Q{n)\ 

2 

=  ylE^  [M((k  +  1)8)2H  -  2M(k8)2H  F  M{(k  -  1)8)2H]  (8.20) 

The  assumptions  apply  certainly  if  Wd  is  a  Levy  subordinator,  i.e.,  an  H'- sssi  1/H'- 
stablc  process,  provided  2 H  <  1  / H' .  Strictly  speaking,  a  subordinator  should  be  a.s. 
increasing  which  is  the  case  for  Levy  motion  provided  H'  >  1.  Note  that  we  do  not 
require  monotonicity  of  Ad.  A  binomial  ‘subordinator’  satisfy  the  assumptions  of  the 
lemma  with  H  —  1/2  since  IEn[Ad&(i)]  =  t.  Finally,  we  note  that  multiplicative  multi¬ 
fractal  processes  with  stationary  increments  have  been  introduced  recently  [70,  74,  105]. 

The  formula  (8.20)  extends  the  usual  auto-correlation  function  of  fGn.  Indeed,  when 
choosing  M.(t)  —  t  it  reduces  to  (7.3),  as  it  should.  In  general,  considering  (8.20),  (7.20) 
and  (7.22),  it  is  probably  fair  to  call 


Hg  = 


TaM(2H)  +  1 
2 


(8.21) 


the  LRD  parameter  of  fractional  Gaussian  noise  in  multifractal  time. 

Consider  the  special  case  H  =  1/2  and  M.  increasing  with  IEn[Ad(t)]  =  t.  Then,  Q  is 
simply  ‘white  (i.i.d.)  Gaussian  noise  in  multifractal  time’.  Its  (2nd  order)  correlations 
vanish  due  to  (8.20).  This  agrees  with  Hg  =  1/2  which  holds  due  to  TM(\)  =  0. 
Nevertheless,  the  process  Q  is  highly  dependent  which  is  revealed  by  the  higher  order 
moments,  i.e.,  by  T&(q)  for  g  /  2.  Mandelbrot  calls  this  the  ‘blind  spot’  of  spectral 
analysis. 
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Simulation 

For  simulations  of  fBm  in  multifractal  time  FB(MF)  it  is  advisable  to  use  a  randomized 
Weierstrass  function 


W(t)  :=  ^  CkXkH  sin(\k M(t)  +  Ak). 

k=  1 

Since  _Ad(s)|  is  for  most  s  and  t  much  smaller  than  \t—s\  (actually  ~  \t— s|“^, 

with  most  a(t)  >  a0  >  1)  any  other  procedure  would  require  high  precision  simulation 
of  fBm  and  might  be  numerically  demanding.  With  C k  i.i.d.  J\f( 0, 1)  and  Ak  i.i.d. 
U  [0,  27r) ,  W  has  zero  mean  increments  with  variance  (compare  [32,  p  248]) 

OO 

JE[(W(t  +  5)  =  2^A2fcHsin2(A k(M(t  +  5)  - 

k= 1 

~  c(M(t  +  <5)  —  M(t))2H. 

This  gives  hope  that  at  least  the  correlation  decay  is  matched. 

Another  simple  case  for  simulation  is  Wiener  motion  in  multifractal  time  WM(MF) 
because  WM,  i.e.,  Bi/2)  has  independent  increments  (compare  Figure  3  (a)-(c)). 

8.4  Pathwise  spectra  for  warped  H- sssi  processes 

In  this  section  we  provide  two  types  of  formulas  for  the  pathwise  spectra  similar  as  be¬ 
fore:  the  ‘multifractal  formalism’  relates  various  spectra  of  subordinated  processes  with 
each  other,  while  the  ‘warp  formulas’  relates  them  to  the  spectra  of  the  subordinator. 

A  first  step  towards  this  goal  is  to  compute  the  deterministic  grained  spectrum  F(a). 
This  is  somewhat  more  involved  than  dealing  with  T(g),  though  we  may  follow  the 
computation  carried  out  for  simple  H- sssi  processes  (compare  (7.11))  with  only  little 
modification.  The  effort  is  worthwhile  since  we  will  gain  information  on  the  pathwise 
spectra  /(a). 

The  strong  monofractal  behavior  of  fBm  relates  the  spectra  of  the  time  change  M. 
and  the  ‘warped’  process  V  in  an  almost  trivial  manner  (see  theorem  8.6)  and  provides 
us  with  a  simple  intuition:  scaling  exponents  multiply  under  subordination.  Since  the 
exponent  of  fBm  is  constant,  i.e.,  h{t)  =  H  for  all  t,  it  is  easy  to  control  the  warped 
spectra.  Allowing  general  H- sssi  processes,  however,  brings  about  the  problem  that 
many  different  combinations  of  singularity  exponents  of  A4  and  V  may  multiply  to 
yield  the  same  exponent  for  V  =  This  will  generally  result  in  ‘smoother’ 

spectra. 

Formulas  and  discussion 

From  the  previous  discussion  we  may  expect  that  the  spectra  of  the  warped  process 
V  and  the  subordinator  Jvl  are  related  through  simple  transforms  which  account  for 
the  aforementioned  ‘smoothing’.  We  first  introduce  and  discuss  these  transforms  on 
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the  spectra  and  establish  then  conditions  under  which  they  equal  the  spectra  of  the 
warped  process.  For  simplicity  we  stick  in  this  section  with  the  case  of  subordinating 
self-similar  processes,  leaving  the  general  processes  for  Section  8.5  when  the  necessary 
formulas  appear  naturally. 

Definition  8.9.  The  H- warped  partition  functions  of  A4  are  defined  as 

(H-sssi:)  ll(*):={I"(9/f)  ^S<4<?’  (8.22) 

and  similarly  for  r]^(g).  Note  that  Tj|4(g)  may  be  — oo  even  for  q  <  q  <q. 

With  the  symbol  H  we  try  to  indicate  the  truncation  of  the  functions.  In  this  notation, 
(8.11)  becomes  simply  Ty(q)  =  T^q),  expressing  that  the  operation  (-)H  gives  the 
correct  passage  from  TM  to  Tv. 

Taking  the  Legendre  transform  of  rjj^,  e.g.,  we  find 

(Tm)*(°)  =  inf  Ua  ~  TmO?))  =  inf  (qa  -  TM(qH ))  =  inf  (V| -  -  TM(q'))  . 

(8.23) 

From  the  last  of  these  expressions  we  learn  that  (rj|/()*(a)  =  r^a/H)  whenever  q  = 
— oo  and  q  =  oo.  This  is  the  case  for  Wiener  motion  Bi/2  when  choosing  the  Holder 
exponents  .  So,  theorem  8.6  tells  us  that  (rj|/1)*(a)  could  indeed  be  a  good  guess 
at  /y  for  general  warped  H- sssi  processes,  but  obviously  only  where  /y  is  concave. 
This  leads  us  to  the  following  definition,  still  in  the  context  of  subordinated  self-similar 
processes  (for  general  processes  see  definition  8.14): 

Definition  8.10.  Assume  that  fM,  respectively  FM,  are  concave.  Then,  the  H- warped 
spectra  of  N4  are  given  by 

(H-sssi,  concave):  f^(a)  :=  (rj[/t)*(a),.  F^(a)  :=  (T|,)*(a). 

(8.24) 

With  the  symbol  A  we  try  to  indicate  the  transformed  shape  of  the  functions  (com¬ 
pare  Figure  9  (b))  as  summarized  in  the  following: 

Lemma  8.11.  Assume  that  fM,  respectively  FM,  are  concave.  Let  fA  and  FA  be  given 
by  (8.24). 

(a)  Shape:  For  a  =  Hr'M(qH )  with  q  <  q  <q 

fM(ya)  =  fM(a/H)  and  (fA)\a)  =  q.  (8.25) 

For  smaller,  resp.  larger  a  the  function  f^(a)  is  linear  with  slopes  q  and  q,  resp. 

(b)  Warping  of  the  Multifractal  formalism:  If  rM(qH)  =  TM(qH )  for  all  q  <  q  <  q, 
then  for  all  a 


(8.26) 
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FIGURE  9.  Left:  To  ‘warp’  a  concave  function  f  with  a  piecewise  linear  function  such  as 
the  grain  spectrum  g  =  Fy  of  a  stable  Levy  motion  one  considers  the  ‘ morphed ’  curves 
9y,z(x)  =  z  —  x(g(y  /  x)  —  1),  singling  out  the  ones  which  touch  f.  Right:  The  resulting  warped 
function  f^(H-)  is  the  set  of  minima  of  those  touching  gy.z.  (Compare  8.4.) 

Proof 

Claim  (a)  follows  from  (8.23).  Claim  (b)  follows  from  (a)  using  corollaries  4.5  and  4.7. 

❖ 

The  following  theorem  establishes  in  what  sense  the  warped  spectra  of  the  subordi- 
nator  A4  equal  the  spectra  of  the  warped  motion.  The  proof  shows  that  the  theorem 
holds  also  for  non-concave  spectra,  thereby  extending  /A  to  the  general  case. 

Theorem  8.12  (Self-similar  processes  in  multifractal  time).  Let  V  be  an  H-sssi 
process  and  AA  an  arbitrary  time  change,  independent  of  V .  Denote  V  :=  V(A4(-)). 
Admit  the  singularity  exponents  o^\  Iifl>  and  v/lf>  under  the  same  conditions  as  in 
theorem  8.5.  Then,  conditioned  on  knowing  AA** 

Fv(a)  \M  =  /£(“)•  (8.27) 

Furthermore, 

Fv(a)  =  F^(a).  (8.28) 

Finally,  if  /m  —  f  M  then  Fy \M  =  Fv \M  (compare  (3.28)  and  (3.26)). 

In  the  special  case  when  V  has  independent  increments  theorem  3.14  allows  to  com¬ 
pute  the  pathwise  grain  spectra  as  well: 

Corollary  8.13  (Levy  stable  motion  in  multifractal  time).  Let  L h  denote  Levy 
stable  motion  and  let  AA  be  an  alm.ost  surely  strictly  increasing  time  change,  indepen¬ 
dent  of  Lh.  Set  V  :=  Lh(A4(-)).  Admit  the  singularity  exponents  cxjf* ,  h ^  and 
under  the  same  conditions  as  in  theorem  8.5.  Assume  that  f  =  /  almost  surely.  Let 

a  be  such  that  f^(a)  >  0  almost  surely.  Then,  for  almost  all  paths  we  have  the  warp 
formula 

fv(a,w )  =  Fy(a)  \M  =  f^(a).  (8.29) 

Moreover,  the  multifractal  formalism  holds  in  the  sense  that 

fv{a,u)  =  Ty(a,  oj)  =  Fv(a)  =  T^(a)  =  (r|,)*(a) 


**In  abuse  of  notation  we  let  Fy(a)  M  =  lim£  limsupn(l/n)  log2  IE[A^(n)(a,  e)  \A4\. 


(8.30) 
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for  any  a  for  which  (a)  ='  /^(a)  >  0. 

Let  us  discuss  some  examples  and  special  cases. 

Example  8.1.  (Binomial  subordinator)  For  a  binomial  time  change  our  results 
are  the  strongest,  since  all  spectra,  i.e.,  all  dimension  and  grain  based  spectra  of  M.h  for 
otjf1 1  tiff'’  as  well  as  w^\  coincide  and  are  strictly  concave  almost  surely.  In  particular, 

f Mb  is  concave,  =  fMb  and  Consequently,  (a)  ='  /^(a)  and 

(8.30)  of  corollary  8.13  applies  to  the  entire  positive  part  of  (a).  4 

Example  8.2.  (Levy  subordinator)  Let  us  assume  that  the  warp  time  is  Levy 
stable  motion  with  Hurst  parameter  H'  while  the  subordinated  process  is  Levy  stable 
motion  with  Hurst  parameter  H ,  i.e.,  V  =  L„(L„,).  Then,  /,  ,  =  /,  ,  =  t!  =  T ! 

- JLj  ti  H' 

by  (7.17)  and  we  may  apply  corollary  8.13.  According  to  (7.13)  these  spectra  are 
piecewise  linear. 

To  compute  fff  ;  let  us  consider  for  simplicity  the  singularity  exponent  efff’ .  In  this 
case  q  =  1/H  and  q  =  —1.  LTsing  that  TLll,(qH )  =  qHH'  —  1  for  — 1  / H  <  q  <  1/(HH') 
— again  by  (7.17)  and  (7.13) —  it  is  straightforward  to  compute  via  (8.24)  and  (8.29) 

(  1  +  (a  —  H  ■  H ')  min(A)  -^L_)  if  a  <  H  •  H’ 

=  {T*  T(a)  =  1  '  if  o  =  H-H'  (8,31) 

y  1  —  (a  —  H  ■  H')  min(l,  A)  if  a  >  H  ■  H' . 

We  will  revisit  this  example  in  the  next  section  and  comment  on  how  to  understand 
this  formula.  4 


Proof  of  theorem  8.12 

As  the  reader  will  notice  the  only  assumptions  on  V  used  in  the  proof  are  stationary 
increments  and  Fy  =  Fv.  We  formulate  the  proof  for  a ^  and  h^f\  Since  only  distri¬ 
butions  matter  here,  (8.12)  is  sufficient  to  deal  with  the  case  of  wavelet  coefficients. 

Let  us  compute  first  the  relevant  probabilities  conditioned  on  M..  For  the  ease 
of  notation  we  assume  that  M.  is  monotonous  so  that  inf  („)  A4(t)  =  M.[k2~n)  and 

sup  (n)  A i(t)  =  M((k  +  l)2“n).  This  simplifies  the  formulas  considerably. 

1k 

We  start  by  looking  at 

Pn(k,  a')  :=  Pn[a  -  e  <  s^,n)(V)  <  a  +  e  \  M,  a'  <  s^\M)  <  a1  +  e']  (8.32) 


Conditioning  on  knowing  M.  we  actually  know  sj^( M.)  (and  not  only  a'  <  s(fl>(A4)  < 
a'  +  e').  Due  to  stationarity  of  increments  of  V  we  have  pn(k,a')  =  pn(0,a').  For 


Sn) 


=  a 


(n) 


we  find  4"  (V)  =  | V(M((k  +  l)2~n)  -  M(k2~n)  )|  =  1^(2“™*  For 

s(F  =  h  'ff  \  s^(V)  becomes  the  supremum  of  increments  of  V  over  an  interval  of  length 


2  nsk  ( M )  adding  only  notational  but  no  logical  difficulties.  Having  expressed  •‘y  (V) 
through  s^\M)  we  may  now  abbreviate  (At)  to  simplify  formulas  without 
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creating  confusion.  A  simple  manipulation  yields 


Pn(k,a')  =  pn(0,a')  = 

_„ew\(a+E)/4’,) 


Pi 


n 


(n)  /  (n)  \ 

<  \V{2~ns>*  )|  <  (wns*  J 


(«)\  (a-e)/Sfcn) 


(8.33) 


a  <  si.  <  a  +  £ 


Using  that  a1  <  <  a'  +  s'  the  asymptotical  behavior  of  pn(k,  a')  is  then 

lirnsup  —  \og2pn(k, a')  =  limsupsj^ — -—log 2pn(k,a')  =  a'  ■  (Fv(a/a')  —  1)  +  o(£,e'). 

n— >oo  n—*oo  77,5!,  ' 

(8.34) 

Here,  o(e,  F)  is  an  error  term  which  goes  to  zero  as  both,  £  and  F  go  to  zero.  If  Fv  =  Fy 
the  liminf„  equals  the  same  expression  (8.34).  We  can  not  conclude,  however,  that  the 
limit  exists  because  pn(k,a')  depends  also  on  e  and  F.  Nevertheless,  we  have  then 
sufficient  control  on  the  asymptotical  behavior  of  pn{k,a!)  for  all  n. 

Let  us  stop  for  a  moment  and  redo  this  calculation  in  the  special  case  of  V  being  an 
H- sssi  process.  We  may  proceed  then  similar  as  in  (7.11)  with  the  only  difference  that 

V  is  now  considered  over  an  interval  of  length  2-nsfc  in  the  case  of  sj^(V)  =  aj^( V ) 
(similarly  for  h^\v)).  Using  stationarity  of  increments  and  self-similarity  of  V  leads 
to 

pn(k,a!)  =  Pn[2~n{a+£)  <\V{2-n^^M))\<2~n{a-£)\  s^\M)} 

=  Pn[n(a  -  H  ■  s^\M)  -  e)  <  s(00\V)  <  n(a  ~  H  ■  s^\M)  +  e)  |  s|,n)(Al)]. 

Using  the  approximation  s^\M)  ~  a'  and  proceeding  as  in  (7.12)  gives 

1  (  (a  —  H  ■  a')q  +  o(s,  F)  if  H  ■  a'  +  e  <  a, 

lim  sup  —  log2  pn{k ,  a')  =  <  o(e,  F)  if  a  —  H  ■  o', 

n^°°  n  y  (a  —  H  ■  a')q  +  o(e,  F)  if  H  ■  a'  —  e  >  a. 

A  moment’s  thought  reveals  that  this  is  indeed  a'  ■  (. Fy(a/a ')  —  1)  +  o(e,F)  as  before. 

Let  us  now  continue  the  proof.  Since  the  asymptotics  of  pn(k,  a ')  does  not  depend 
on  k  we  may  write  pn(lF)  =  pn(k,l£')  and  group  terms  as  follows 

2n  — 1 

-Jin  :=  ^2  -  £  <  4n)(V)  <  a  +  e  I  M\  =  y^Pn^le')  ■  N£i(lF,F/2).  (8.35) 

k= 0  l 

This  approach  closely  follows  the  techniques  of  the  proof  of  theorem  4.2  and  lemma  3.6. 
First,  one  notices  that  for  fixed  e  and  e'  there  are  only  finitely  many  indices  l  to  consider. 
Hence  the  asymptotic  behavior  of  (8.35)  is  dominated  by  the  asymptotically  largest 
term  appearing  in  the  sum. 

The  asymptotics  of  pn  (8.34)  combined  with  the  definition  of  fM  (1.8)  yield 
lirnsup  -  log2  [pn(a')  ■  N^A(a',  F /2)]  =  /»(a')  +  a'  ■  (Fv(a/a')  -  1)  +  o(e,  F).  (8.36) 

n— »oo  ^ 
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As  in  (3.16)  and  (4.4)  one  concludes  that 
1 

lim  sup  —  log2  Tin  =  sup  ( /x(a ')  +  a  ■  ( Fy(a/a ')  —  1))  +  o(e,  e'). 

n—>oc  Tl  a' 

For  this  computation  one  uses  the  general  fact  that  lim  inf(xn)+lim  inf  (yn)  <  lim  inf  (xn+ 
yn)  <  lim  inf  (xn)  +  limsup(yn)  <  limsup(xn  +  yn)  <  limsup(xn) +  limsup(yn).  Because 
some  of  these  inequalities  may  be  strict,  however,  one  needs  to  know  that  at  least  one 
°f  fM  =  lM  or  Fv  =  Fv  holds  to  obtain  this  limsup;  and  both  equalities  have  to 
hold  if  one  likes  to  estimate  the  liminf.  For  H- sssi  processes,  Fy  =  F_v  follows  from 
corollary  7.5.  Letting  £  and  s'  tend  to  0  we  get,  finally, 

Fv(a)  I  m  =  sup  (/m(o')  +  a!  •  ( Fv{a/a! )  -  1))  (8.37) 

a' 

Formula  (8.37)  captures  the  multifractal  interaction  between  warp  time  and  warped 
motion.  It  makes  most  explicit  how  the  combination  of  singularity  exponents  of  A4 
and  V  results  in  the  spectrum  of  V  =  V(M).  The  RHS  of  (8.37)  is  the  announced 
definition  of  /f)  in  all  generality  (compare  (8.24)).  We  will  discuss  its  properties  as 
well  as  the  fact  why  we  consider  it  a  transform  of  /m  rather  than  of  Fv  momentarily. 

The  proof  of  the  first  part  of  theorem  8.12  is  complete  now.  For  the  second  part,  we 
write  by  the  law  of  total  probability,  by  (3.22)  and  by  rearranging  terms, 

2n— 1  2n  — 1 

£  C„[  |4”>(V)  -  a\  <  e]  =  £][>„(fe')  ■Pa[\sf>(M)-le'\<e'l2\ 

k= 0  k= 0  l 

=  £p„(i£')'En[Jvy(k',e72)] 

l 

Comparing  with  (8.35)  and  recalling  (3.23)  it  is  clear  that  the  same  analysis  as  above 
applies,  with  /m  replaced  by  F,vf .  This  proofs  (8.28).  <) 

Proof  of  corollary  8.13 

To  proof  the  corollary  we  apply  theorem  3.14  conditioned  on  knowing  M..  Note  that  the 
singularity  exponents  s[n^  =  \LH(M({k  +  l)2-n))  —  LH(M(k2~n))\  (k  =  1,2, . . .)  be¬ 
come  independent  when  conditioned  on  A4.  Similarly  for  s[n>  =  sup{|L^(u)  —  LH(u')  : 
M((k  —  1)2-")  <  u  <  u'  <  M((k  +  2)2~n)|,  at  least  in  the  sense  of  remark  3.15. 
Finally,  (8.30)  follows  from  the  multifractal  formalism  corollary  4.1  applied  to  V  and 
from  fv  =  f£l  =  Tlli  =  Tv.  0 

General  definition  of  /A 

Motivated  by  formula  (8.37)  we  extend  definition  8.10: 

Definition  8.14.  The  F-warped  spectra  of  JGl  are  given  by 

f^(a)  =  sup  {/mW)  +  a'  ■  ( Fy(a/a ')  -  1)) 


(8.38) 


and  similar  for  F ^ 
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Let  a  mark  a  value  for  which  FY  attains  its  maximal  value,  which  is  always  1.  Then, 
(FY(a/a!)  —  1)  =  0  at  a'  =  a/a  which  implies  that 

/£(<*)  >  fM (?)• 

CL 

In  order  to  get  a  more  graspablc  idea  of  fA  consider  the  two-parameter  family  of 
functions  ga,z(a')  ■  =  z  —  a1  ■  (. FY{a/a ')  —  1).  First,  observe  that  each  such  function  has 
as  minimum  point  in  (a/a,z).  Second,  since  /x(a')  <  1  we  have  ga>z  >  Jm  for  z  >  1. 

We  are  now  ready  to  give  an  explanation  of  (8.38)  as  a  ‘touching’  or  ‘morphological’ 
transform  (compare  Figure  9): 

/£(<*)  =  min{2  :  9a, . z  >  Im}-  (8-39) 

In  other  words,  z  =  f%(a)  is  the  critical  value  of  z  for  which  the  curves  g(l)Z  and  Jm 
touch. 


Example  8.3.  (Warped  fBm)  In  the  special  case  of  Y  =  F>h  and  choosing  the 
Holder  exponents  /i[n) ,  the  spectrum  FY  reduces  to  the  point  (a,  1)  and  we  hnd  as 
earlier 

fu(a)  =  fM(a/a )  (if  q  =  ~oo,  q  =  oo.)  (8.40) 

Example  8.4.  (Warped  H- sssi  process)  For  an  //-sssi  process  Y  and  general  s(/L) 
one  has  a  =  H  and  by  (7.13) 


J  z  —  (a  —  H  ■  a!)q  for  a'  <  a, 
\  z  —  (a  —  H  ■  a')q  for  a'  >  a. 


(8.41) 


Note  that  q  <  0  <  q.  Thus,  gatZ  has  indeed  the  shape  of  Figure  9  (a),  with  a  unique 
minimum  in  (a/ H,  z)  and  with  slopes  F[q  and  Hq  to  the  left  and  right  of  this  minimum, 
respectively. 

As  we  have  seen  in  (8.39)  one  needs  to  draw  the  ‘touching’  gaz  and  to  follow  their 
minima.  Keeping  the  simple  form  of  ga  z  in  mind  we  may  conclude  as  follows: 

If  a  is  such  that  the  appropriate  gajZ  is  touching  fj ^  with  its  minimum,  then  we  have 
/m(«)  =  fMia/H). 

For  all  other  a  the  appropriate  ga^z  is  touching  /m  with  one  of  its  linear  parts,  the 
small  (large)  a  on  the  one  of  slope  Hq  (Hq).  Using  (8.39)  and  z  =  ga>z(a/H)  it  follows 
that  //)(a)  =  z  depends  linearly  on  a,  for  small  (large)  a  with  slope  q  (q). 

Intuitive  recipe  Since  an  H- warped  f'M  is  continuous,  computing  it  amounts  to 
finding  ‘tangents’  to  /m  of  slope  Hq  and  Hq  on  the  left  and  right,  respectively.  Fol¬ 
lowing  the  graph  of  /m  in  its  center  part  and  the  tangents  at  the  outskirts  on  gets  a 
plot  of  /£(#■). 

Rigorous  recipe  One  takes  the  Legendre  transform  of  fj^(-/H)  (which  is  nothing 
but  tm('H)),  cuts  off  the  slopes  q  outside  the  interval  (q,q)  and  transforms  back.  Thus, 
(8.38)  is  a  true  extension  of  (8.24).  4k 
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Note  that  /A  is  not  always  concave,  though  it  may  remove  some  of  the  non-concave 
parts  of  /.  For  toy  models  of  non-concave  /  which  will  lead  to  non-concave  /A  see 
[91,  94], 

Finally,  let  us  explain  why  (8.38)  should  be  viewed  as  a  transform  of  / m  rather  than 
one  of  Fy ■  First,  the  simplicity  of  (8.41)  and  (8.24)  may  be  misleading.  In  general, 
the  smoothing  ‘kernel’  a'  ■  ( Fy(a/a ')  —  1)  will  change  its  shape  as  a  changes.  A  second 
reason  is  found  in  (8.40).  Third,  Fy  is  deterministic  and  f m  is  not. 

Example  8.2.  (Levy  subordinator,  revisited)  Let  us  apply  formula  (8.38)  to  the 
situation  introduced  in  Example  8.2  on  page  76  and  compute  the  positive  part  of 
.  Recall  that  q  —  l /Ft  and  q  —  —  1  for  the  singularity  exponent  a ^  considered 
there.  According  to  (7.17)  the  graph  of  f^  (-/H )  consists  of  two  line  segments  joining 
(HH1, 1)  with  the  origin  and  with  (77(77'  +  1) , 0)  (slopes  1/(7777')  and  —1/H). 

At  this  point  it  is  already  clear  that  fy(HH')  =  (7777')  =  fLH,(H')  =  1  which 

marks  the  peak  of  fy.  This  does  not  come  as  a  surprise:  the  (Lebesgue)  almost  every¬ 
where  present  Holder  exponents  F[  and  77'  of  motion  and  time  change  (see  [54,  Rem.  4]) 
combine  into  the  (Lebesgue)  almost  everywhere  present  Hdlder  exponent  7777'  of  the 
compound  process. 

Let  us  assume  77'  >  1  in  order  to  have  an  increasing  warp  time.  Then,  fy  =  f£  = 
fLH,(-/H)  in  their  left  part,  i.e.,  for  a  <  HH' .  This  is  because  their  slope  is  there  equal 
to  1  /{HH')  which  lies  in  the  interval  (g,  q)  =  (—1,1  /H).  Consequently,  fy  will  pass 
through  the  origin. 

If  H  >  1,  then  fy  =  f£  =  hH,(-/H)  also  in  the  right  part  since  the  slope  is  there 
—  1  /H  which  is  within  (g,g)  =  (—1,1/77): 


V  —  LH>l(LHI>1)  :  /v(a)  —  j  ^  _ 

However,  if  77  <  1  then  the  right  part  of  f£ 
allowed: 

V  =  Lff<i(Lff/>i)  :  fv(a )  =  j  ^  ^  |a 


77  ■  77')/(77  ■  77') 
77- 770/77 

will  have  slope 


77  ■  T7')/(T7  ■  H') 
77-T7') 


if  0  <  a  <  77  ■  H' 
if  a  >  77  ■  H'. 

1,  the  smallest  one 


if  0  <  a  <  77  ■  H' 
if  a  >  77  ■  77'. 


With  (8.38)  we  have  now  a  better  understanding  of  the  difference  of  these  two  cases. 
First,  when  77  <  1  the  subordinated  process  Lh{Lh')  is  somewhat  more  regular  than 
a  comparison  with  ffim  would  make  one  think,  i.e.,  f£  >  fLH,{-/H)  for  a  >  HH' 
(note  that  this  agrees  with  (8.25)).  For  an  explanation  note  that  the  ‘touching  point’ 
of  fjJH,  with  a\FL[1(a/a')  —  1)  occurs  at  the  peak  of  fL  meaning  that  the  (Lebesgue) 
almost  everywhere  present  L///-points  with  h{t)  =  H'  >  1  control  the  entire  right  side 
of  the  multifractal  spectrum  of  V.  The  left  part  is  governed  by  the  (Lebesgue)  almost 
everywhere  present  Lj/-points  with  h{t)  =  77.  Second,  in  the  case  77  >  1  the  entire 
spectrum  is  controlled  by  the  same  (Lebesgue)  almost  everywhere  present  L//-points 
with  h(t)  =  77.  4|h 
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8.5  Multifractals  in  Multifractal  Time 

In  Section  8.4  we  actually  established  theorem  8.12  not  only  for  H- sssi  processes  V 
but  for  general  processes  Y  with  stationary  increments  and  some  regular  multifractal 
scaling: 

Theorem  8.15  (Multifractals  in  multifractal  time).  LetY  and  M.  be  independent 
processes  and  set  y(t)  =  Y(A4(t)).  Assume  that  the  distributions  ofY(t  +  5)  —  Y(t) 
are  independent  of  t.  Admit  the  singularity  exponents  oi^\  h!f'}  and  under  the 
same  conditions  as  in  theorem  8.5.  Assume  that  Fy  =  FY.  Then,  we  have  the  warp 
formula 

Fy(a)\M  =  fM(a)  and  Fy(a)  =  F^(a),  (8.42) 

and  A 

Ty{q)  \M  —  TM{Ty{q)  +  1)  Ty(q)  =  TM(TY(q)  +  1).  (8.43) 

Moreover,  if  fM  =  f_M  then  Fy \M  =  Fy \M. 

Proof 

To  obtain  Ty  we  simply  take  the  Legendre  transform  in  (8.42)  .  Assuming  a'  >  0  we 
find 

Ty(q)  m  =  Fy(q)  M  =  inf  (qa  -  Fy (a)  =  inf  (qa  -  /£(a)) 

=  inf  (qa  +  inf  (— /yAa')  —  a'lFyiala1)  —  1))) 

=  'sf  + (* + “  (r6  -  Fy0))) 

=  inf  (-/At(a')  +  a'  (1  +  Ty(q ))) 

a' 

=  rMa+Ty(i))  =  TMa+Ty(q)). 

With  similar  manipulations  one  obtains  Ty.  As  a  matter  of  fact,  (8.43)  can  easily  be 
obtained  directly  using  the  technique  of  (8.33).  <0> 

Formula  (8.43)  shows  how  to  extend  definition  8.9: 

Definition  8.16.  The  Y -warped  partition  functions  of  M.  are  given  by 

Tl(q )  :=TM(TY(q)  +  l)  (8.44) 

and  similarly  for  (q) . 

d  Formula  (8.43)  reads  the  same  for  any  c-ary  — not  only  dyadic —  scaling  (compare  foot¬ 
note  t).  If  we  replaced,  however,  S^n\q)  by  IEn  2~nqsk  ,  =  c~dnS^n\q)  in  the  definition  of 
T(q)  (compare  (3.11)  and  (3.19))  then  (8.43)  would  become  the  appealing 

TY{M)  =  Tm(Ty). 
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Indeed,  since  for  H- sssi  processes  V  we  have  Ty(q)  =  qH  —  1  for  q  <  q  <  q  and  may 
indeed  write  (8.22)  in  the  form  T_j|/1(g)  =  TM{Tv(q)  +  1)  in  that  case. 

Example  8.5.  (Subordinated  Multifractals)  Turning  around  the  usual  roles  we 
may  consider  y  =  A4(LH)  and  get 


TM(LH)(q ) 


H  ■  Tm  (q)  +  H  —  1  if  q  <  TM  (q)  +  1  <  q 
—oo  else. 


Assuming  differentiability  this  leads  with  (4.11)  to  Fm(lh){o)  =  HFm(ci/H)  +  1  —  H 
which  is  valid  for  a  =  HT'M{q)  with  q  in  the  above  range;  for  other  a,  FM(Lh)(cl)  will 
be  linear  of  slope  such  as  to  make  the  graph  smooth  for  all  a.  4 

Combining  theorem  8.15  and  3.14  we  get: 

Corollary  8.17  (Levy  processes  in  multifractal  time).  Let  Y 

theorem  8.15  with  the  additional  assumption  that  the  increments  ofY 
Then,  we  have  the  pathwise  warp  formula 

fy(a,v)  \M  =  Fy(a)  \M  =  f*(a).  (8.45) 

Moreover,  the  multifractal  formalism 

fy{a,u)  =  Ty(a,uj)  =  T$(a)  (8.46) 

holds  for  any  a  for  which  (a)  ='  /^(a)  >  0. 

A  process  Y  with  stationary  and  independent  increments  as  in  corollary  8.17  is 
called  a  Levy  process.  In  the  special  case  when  Y  is  7-stable,  it  is  necessarily  l/7-sssi, 
in  other  words  Levy  stable  motion,  which  we  dealt  with  earlier.  Jaffard  in  [54]  has 
obtained  the  almost  sure  dimension  based  spectrum  for  a  large  class  of  Levy  processes, 
in  particular  Levy  stable  motion.  Here,  we  computed  the  grain  based  spectra  of  warped 
Levy  processes,  provided  one  has  knowledge  on  the  moments  or  the  tail  behavior  of 
the  marginals  of  Y  in  order  to  compute  Ty  or  Fy. 

Example  8.6.  (Levy  motion  in  multifractal  time)  A  special  case,  in  which  all 

assumptions  of  corollary  8.17  hold,  is  present  if  the  warped  process  is  Levy  stable 
motion,  i.e.,  Y  =  Ljj,  and  where  the  warp  time  is  a  binomial  cascade,  i.e.,  At  =  A4b- 
This  is  true  for  the  singularity  exponents  h)”'1  and  a\ without  further  ado,  while  for 
w one  has  to  check  the  existence  of  moments  of  negative  order  for  An)k  (compare 
corollary  6.2).  4k 


and  M.  be  as  in 
are  independent. 


Dimension  based  spectra 

There  is  a  simple  class  of  multifractals  in  multifractal  time  on  [0, 1]  for  which  the  full 
multifractal  formalism,  in  particular  the  dimension  based  spectra  are  known.  As  we 
will  see  the  same  warp  formula  (8.43)  will  hold  despite  the  fact  that  strictly  speaking 
the  example  does  not  fit  in  the  above  framework. 
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Full  multifractal  warp  formalism  Let  Mb  and  /Q>  be  independent  binomial  dis¬ 
tribution  functions  as  in  (5.2)  with  properties  (i)-(iii)  (see  (5.5))  and  with  multipliers 
M^>}  and  K[n\  respectively.  To  ensure  stationarity  at  least  of  their  dyadic  increments 

assume  in  addition  that  M0  =  Mi  and  K0  =  K\ . 

Now,  let  K-l  be  the  inverse  function  of  /Cj,  just  as  in  Section  5.6,  i.e. ,  Kb{KL\{t))  =  t. 
Let  y  =  Mb(ICl).  Since  Mb  and  1C\  are  both  continuous  and  strictly  increasing,  so 
is  y  and  we  may  define  a  measure  by  setting  /x([0,  t])  :=  y(t).  Let  us  check  at  how 
p  distributes  mass  on  the  unit  interval  [0,1].  Therefore,  let  tn ^  :=  Kb[k2~n).  These 
points  partition  the  unit  interval  into  2"  intervals  :=  [tn,k,tn,k+ 1]  which  have  the 
length  ■/)'“  =  tn<k+1  -  tn<k  equal  to 

ICMK  +  1)2-”)  -  !Cb(kn 2-")  =  At>  ■  At:,11  •  •  •  A'W  •  A<°> 

and  /(-mass  =  Y{tn,k+i)  —  Y'(tnk)  equal  to 

MMK  +  1)2-“)  -  Mb(K 2-“)  =  «Sf  •  M'”-1’  ■  ■  ■  M<;> .  Mq°\ 

Thus,  we  recognize  in  this  measure  /x  immediately  the  measure  /ic  of  Section  5.6  with 
c  =  2  and  Li'1'1  =  K\n\  In  other  words,  p  is  a  special  case  of  a  so-called  statistically 
self-similar  measure  [33,  6,  8]  for  which  equality  in  the  multifractal  formalism  for  the 
singularity  exponents  cqn'1  and  /i[n)  is  well  established,  and  for  which  T  is  given  through 
the  implicit  equation  (compare  (5.41)) 

E  \M%K~T{q)  +  M?K~T{q) ]  =  2E  [. MqK~T{q) ]  =  2E  [Mq]  E  [ii  "r(<7)]  =  1.  (8.47) 

For  a  comparison  of  (8.47)  with  the  usual  warp  formula  (8.43)  let  us  recall  that 

E[M9]  =  2~TMb^-1  and  E[A'^(9 }]  =  29'-1  (8.48) 

from  (5.25)  and  (5.40).  Setting  q1  =  T_Aa6(g)  + 1  in  the  last  expression  we  see  that  (8.47) 
is  indeed  solved  by  T(q )  =  T/ct(T_A46(g)  +  1).  We  see  a  reason  for  (8.43)  to  hold  in  the 
fact  that,  though  Mb  does  not  have  stationary  increments  of  all  lengths,  nevertheless 
the  lengths  |J^|  and  the  masses  ^(J^)  are  stationary,  i.e.,  of  equal  distribution  for 
k  =  0, . . . ,  2n  —  1  (n  fixed). 

Remark  8.18.  (Independence  of  warp  time  and  warped  motion)  Note  that 
(8.43)  does  not  necessarily  hold  if  Y  and  M  are  dependent.  This  follows,  e.g.,  when 
choosing  /Cf,  to  be  pathwise  equal  to  Mb-  Then,  we  get  y(t)  =  t  and  Ty(q)  =  q  — 
1.  Combining  Tjcb(—TKt(q))  =  —  q  (which  follows  from  (8.48))  with  (8.43)  we  find 
Ticb( 1  ~  q)  =  TKb(-Ty(q ))  =  ~TMb(q)  -  1  =  - TKb(q )  -  1.  Setting  q  =  1/2  gives 
T/cb(  1/2)  =  —1/2.  Since  7)c6(0)  =  —1  and  Tfch ( 1 )  =  0  we  conclude  that  T/cb  must  be 
linear  for  0  <  q  <  1  due  to  concavity.  But  this  does  not  hold  in  general  and  we  must 
conclude  that  (8.43)  may  fail  when  the  components  are  dependent.  & 
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Warped,  inverse  and  relative  multifractals 


In  the  framework  of  so-called  relative  multifractal  analysis  one  aims  at  obtaining  more 
flexible  multifractal  spectra  by  replacing  the  usual  singularity  exponents  of  a  process 
A4,  say 


W(M) 


log  | M((k  +  1)2“")  -  M(k2~n)\ 

!og|4n)l 


log  \M((k  +  1)2“")  -  M(k2~n)\ 
—n  log(2) 


by  exponents  which  take  into  account  a  different,  hopefully  more  sensitive  measure  of 
I  If1  than  its  length.  More  explicitly,  the  multifractal  spectra  of  M.b  relative  to  JCb  are 
in  terms  of  the  singularity  exponents 


^  (n)  =  log  \Mb((k  +  l)2~n)  —  Mb[k2~n)\ 
fk  log  \1Cb((k  +  l)2_n)  —  JCb(k2~n)\ 


In  the  notation  of  [95]  rMb^K'b(9)  =  limn  rjC  (kn 2  n  — >  6).  Comparing  with  y  = 
M.b(K\)  from  above  we  find  immediately  that  cy^\y)  =  log  /i(jjf^)/  log  coincides 
with  the  ‘relative  exponent’  r\f\M.b/ ICb) .  The  only  difference  is  that  they  ‘live’  in 
two  different  spaces:  rbl>  is  defined  for  ^-intervals  lj^\  and  for  f-intervals  Jbl)  = 


Nevertheless,  we  may  immediately  say  that  the  number  of  intervals  of  length 
2~n  with  r["')  ~  a  grows  as  (2~n)~fla'>  (compare  (3.5))  where  /  =  fMb(tf)  =  T*  from 
(8.47).  This  is  exactly  the  formula  ‘/r(o:)  =  /?*(«)’  of  [95,  p.  160]  and  can  also  be  found 
in  [63]. 

However,  it  might  be  appropriate  to  abandon  the  Euclidean  length  ijf1'*  in  0-space 
altogether  and  to  replace  it  by  \Kb{{k  +  l)2“n)  —  K,b[k2~n)  \  not  only  what  concerns  the 
singularity  exponent  but  also  in  the  spectra.  Again,  (8.37)  is  key.  It  helps  understand 
how  to  do  so  as  it  shows  how  to  transfer  resolution  from  one  space  to  the  other:  the 
value  of  a'  for  which  the  supremum  is  attained  is  the  predominant  distortion  factor  by 
which  K)h  ‘expands’  dimension.  This  suggests  the  formula  /)»  Mb(a)  =  a!  ■  fMb^)(a)- 

It  holds  indeed  in  the  form  ‘dim (K^ib^ICb)  =  /3*(a )  •  aKb(q)'  of  [95,  p.  160]. 

Choosing  ICb  to  be  iid  with  M.b  this  formula  suggests  what  to  expect  of  a  multifractal 
analysis  of  a  random  measure  relative  to  an  independent  realization  of  itself. 

Understanding  the  multifractal  structure  of  M.b  with  respect  to  lCb  is,  thus,  essentially 
equivalent  to  understanding  the  structure  of  M.b  warped  by  the  inverse  multifractal 
of  fCb.  Since  inverse  measures  can  be  defined  in  general  [94]  and  since  the  basis  of  the 
relative  multifractal  formalism  has  been  expanded  [17]  this  principle  is  bound  to  hold 
in  far  broader  context. 


Conclusion 

We  have  presented  the  coarse  grain  aspects  of  multifractal  analysis  from  a  most  general 
point  of  view  which  allows  to  connect  the  various  different  approaches  currently  ern- 
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ployed  in  the  literature,  i.e.,  various  singularity  exponents  as  well  as  different  notions 
of  dimensions  and  scaling. 

The  framework  presented  proves  convenient  for  exploring  the  multifractal  properties 
of  the  well  known  random  cascade  measures  as  well  as  self-similar  processes.  But  it 
applies  also  to  the  novel  compound  processes  we  introduced.  In  subordinating  fractional 
Brownian  motion  (ffim)  -  and  more  generally  stable  motion  -  to  multifractal  time  we 
have  only  scratched  the  surface  of  a  new  class  of  multifractal  objects.  Many  questions 
must  jump  to  mind  immediately  which  are  beyond  the  scope  of  this  paper. 

A  problem  of  general  interest  is  the  construction  of  multifractal  measures  with  sta¬ 
tionary  increments.  In  [14,  70,  74,  105]  different  classes  of  such  processes,  each  with 
appealing  properties,  are  proposed. 

From  a  multifractal  point  of  view  one  can  ask  for  more  precise  spectral  information 
such  as  exact  Hausdorff  gauge  functions  and  Hausdorff  measures  of  sample  paths, 
an  issue  which  has  been  solved  for  graphs  and  images  of  fBm  only  in  the  last  years 
[107,  108], 

In  view  of  the  dimension  based  multifractal  spectra  it  would  be  interesting  to  know 
to  what  extent  equality  in  the  multifractal  formalism  holds.  A  first  approach  to  this 
problem  could  be  to  focus  on  time  change  of  Levy  processes. 

From  a  simulation  or  engineering  point  of  view  the  new  models  raise  the  question 
of  efficient  and  accurate  estimation  of  parameters  and  synthesis.  There  is  a  general 
need  for  reliable  empirical  evidence  of  multifractal  properties  of  time  series,  as  well 
as  for  plausible,  physical  explanations  of  multiplicative  structures  towards  better  un¬ 
derstanding.  Research  in  this  area  could  profit  from  the  large  literature  discussing  the 
extraction  of  multifractal  parameters  out  of  physical  systems  which  took  its  beginning 
with  [42,  45,  46].  On  the  statistics  of  multifractal  estimators  from  time  series  some 
work  has  been  done  recently  [41,  3,  40]. 

The  modelling  of  fractal  and  multifractal  processes  has  so  far  relied  mostly  on  fBm 
and  cascades.  For  TCP  traffic  loads  and  web  access  hie  size  processes,  e.g.,  multifractal 
cascades  such  as  the  binomial  provide  accurate,  yet  parsimonious  models  and  have  had 
a  considerable  impact  [89,  96,  38,  88].  Other  traffic  loads  which  are  not  subject  to  large 
time  delay,  such  as  video,  are  apparently  well  modelled  by  fBm.  Finally,  for  financial 
data,  fBm  in  multifractal  time  has  been  proposed  as  a  model  in  [68]. 

For  simulation  purposes  of  fBm  in  multifractal  time  a  good  idea  is  to  start  with 
Wiener  motion  which  has  independent  increments.  Alternatively,  one  may  submit  the 
randomized  Weierstrass  function  to  multifractal  time  since  it  is  more  simple  to  synthe¬ 
size  than  fBm.  Moreover,  modifying  approximate  models  based  in  the  wavelet  domain 
such  as  WIG  (see  [88])  or  the  FFT  based  models  of  [20]  should  be  more  practical  than 
direct  synthesis  in  the  time  domain  in  the  presence  of  the  strong  correlations  of  fBm. 
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